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ELLIPTIC  AND 
PARABOLIC  PARTIAL 
DIFFERENTIAL 
EQUATIONS  IN 
NUMERICAL  GRID 
GENERATION 


1  INTRODUCTION 


The  primary  objective  of  this  study  is  to  explore  the  feasibility  of  using 
parabolic  partial  differential  equation  techniques  for  numerical  grid  genera¬ 
tion  for  two-dimensional  aerodynamic  configurations.  The  contents  include 
a  discussion  of  grid  generation  concepts  and  schemes  in  the  literature,  iter¬ 
ative  methods  for  numerical  grid  generation  and  two  differential  grid  gen¬ 
eration  schemes:  (1)  an  elliptic,  and  (2)  a  parabolic  scheme.  A  detailed 
mathematical  and  numerical  representation  of  both  schemes  is  given.  The 
main  purpose  in  the  treatment  of  the  elliptic  scheme  is  to  introduce  the 
necessary  transformations  to  provide  equations  that  will  establish  the  par¬ 
ticular  parabolic  scheme  development. 

The  grid  gpneration  scheme  is  derived  from  a  pair  of  model  parabolic 
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equations  in  which  the  coefficients  are  determined  by  use  of  a  related  set 
of  two-  dimensional  elliptic  partial  differential  equations.  Coefficients  that 
have  effects  similar  to  the  exponential  forcing  functions  for  grid  spacing 
control  in  elliptic  grid  generation  are  determined.  This  formulation  involves 
discretization  of  the  governing  equations  on  a  quasi-uniform  grid  in  the 
computational  plane  and  using  linear  interpolation  between  the  current 
grid  line  and  outer  boundary. 

The  basic  idea  in  grid  generation  by  use  of  parabolic  PDE’s  is  to  specify 
the  body  surface  as  an  “initial”  boundary,  and  treat  the  outer  boundary 
as  a  constraint.  Then  coordinate  lines  are  generated  by  marching  outward 
from  the  body  surface. 

The  parabolic  grid  generation  method  is  extended  to  multicomponent 
airfoils,  arbitrary  multiple  bodies  and  cascade  grid  generation  for  the  first 
time.  We  present  representative  results  as  a  measure  of  the  efficiency  of  the 
parabolic  scheme. 

2  GRID  GENERATION 

Numerical  grid  generation  has  become  an  integral  part  of  computational 
fluid  dynamics  (CFD)  and  is  one  of  the  most  important  topics  in  the  de¬ 
velopment  of  flow  solutions  in  complex  flow  domains.  The  grid  here  is 
understood  to  be  an  organized  set  of  points  formed  by  the  intersections 
of  the  coordinate  lines  of  a  numerically  generated  boundary-conforming 
curvilinear  coordinate  system  having  the  same  dimensionality  as  the  phys¬ 
ical  region.  The  main  feature  of  such  a  system  is  that  some  coordinate 
line  (surface  in  three  dimensions)  is  coincident  with  each  segment  of  the 
physical  boundary.  There  are  basically  two  decision  stages  involved  in  the 
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Part  I 

ELLIPTIC  AND 
PARABOLIC  PARTIAL 
DIFFERENTIAL 
EQUATIONS  IN 
NUMERICAL  GRID 
GENERATION 


1  INTRODUCTION 

■'The  primary  objective  of  this  ^tudy  is  to  explore  the  feasibility  of  using 
parabolic  partial  differential  equation  techniques  for  numerical  grid  genera¬ 
tion  for  two-dimensional  aerodynamic  configurations.  The  contents  include 
a  discussion  of  grid  generation  concepts  and  schemes  in  the  literature,  iter¬ 
ative  methods  for  numerical  grid  generation  and  two  differential  grid  gen¬ 
eration  schemes:  (1)  an  elliptic,  and  (2)  a  parabolic  scheme.  A  detailed 
mathematical  and  numerical  representation  of  both  schemes  is  given.  The 
main  purpose  in  the  treatment  of  the  elliptic  scheme  is  to  introduce  the 
necessary  transformations  to  provide  equations  that  will  establish  the  par¬ 
ticular  parabolic  scheme  development.  ' '  .  ' 

The  grid  generation  scheme  is  derived  from  a  pair  of  model  parabolic 
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equations  in  which  the  coefficients  are  determined  by  use  of  a  related  set 
of  two-  dimensional  elliptic  partial  differential  equations.  Coefficients  that 
have  effects  similar  to  the  exponential  forcing  functions  for  grid  spacing 
control  in  elliptic  grid  generation  are  determined.  This  formulation  involves 
discretization  of  the  governing  equations  on  a  quasi-uniform  grid  in  the 
computational  plane  and  using  linear  interpolation  between  the  current 
grid  line  and  outer  boundary. 

The  basic  idea  in  grid  generation  by  use  of  parabolic  PDE’s  is  to  specify 
the  body  surface  as  an  “initial”  boundary,  and  treat  the  outer  boundary 
as  a  constraint.  Then  coordinate  lines  are  generated  by  marching  outward 
from  the  body  surface. 

The  parabolic  grid  generation  method  is  extended  to  multicomponent 
airfoils,  arbitrary  multiple  bodies  and  cascade  grid  generation  for  the  first 
time.  We  present  representative  results  as  a  measure  of  the  efficiency  of  the 
parabolic  scheme. 

2  GRID  GENERATION 

Numerical  grid  generation  has  become  an  integral  part  of  computational 
fluid  dynamics  (CFD)  and  is  one  of  the  most  important  topics  in  the  de¬ 
velopment  of  flow  solutions  in  complex  flow  domains.  The  grid  here  is 
understood  to  be  an  organized  set  of  points  formed  by  the  intersections 
of  the  coordinate  lines  of  a  numerically  generated  boundary-conforming 
curvilinear  coordinate  system  having  the  same  dimensionality  as  the  phys¬ 
ical  region.  The  main  feature  of  such  a  system  is  that  some  coordinate 
line  (surface  in  three  dimensions)  is  coincident  with  each  segment  of  the 
physical  boundary.  There  are  basically  two  decision  stages  involved  in  the 
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discretization  of  a  flow  field.  The  first  involves  a  decision  about  the  grid- 
generation  concept  to  be  used,  and  the  second  involves  a  decision  about  the 
grid  generation  scheme  to  be  employed.  Several  grid-generation  concepts 
and  schemes  are  summarized  in  fig  1. 

Most  discretization  methods  can  be  classified  as  based  upon  one  or  more 
of  the  concepts  in  fig.  1.  These  include  single-module,  multi-block,  interfer¬ 
ing,  component-adaptive  overlapping,  and  component  adaptive  interfacing 
approaches.  Probably  the  easiest  and  most  popular  concept  is  the  single 
module,  in  which  the  discretized  flow  domain  is  transformed  into  a  single 
computational  rectangle  (in  2-D)  or  cube  (in  3-D). 

The  multi-block  approach  corresponds  to  linking  together  several  su¬ 
ch  blocks.  Two  methods  are  discussed  by  Rubbert  and  Lee  [1].  The  first 
method  is  to  generate  the  grid  separately  within  rectangular  subdomains 
(blocks).  Some  of  the  block  boundary  surfaces  no  longer  correspond  to 
boundary  surfaces  of  the  original  problem,  but  instead  separate  adjacent 
blocks.  These  are  called  field  boundaries.  Solution  of  the  grid  genera¬ 
tion  equations  requires  grid  boundary  conditions  on  these  field  boundaries. 
Grids  generated  in  this  manner  are  termed  “patched  grids”  or  “patched 
coordinate  systems”.  The  second  method  is  to  solve  the  grid  generation 
equations  in  the  entire  block-  structured  domain  as  a  single  grid  genera¬ 
tion  problem.  In  this  case  grids  will  be  analytical  across  field  faces,  but 
solution  of  the  problem  is  more  complex.  Grids  generated  in  this  manner 
are  called  “directly  solved  multi-  block  grids”.  An  example  of  an  effective 
multi-block  approach  is  given  by  Thomas  [2]  for  discretizing  a  wing  fuselage 
configuration. 

Once  the  desired  grid-generation  concept  has  been  determined  for  a 
given  flow  problem,  the  next  step  is  to  determine  the  actual  location  of 
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the  grid  points.  To  do  this,  a  variety  of  grid  generation  schemes  are  avail¬ 
able.  Basically,  the  grid-generation  schemes  for  curvilinear  grids  are  of  two 
general  types  (fig.  2)  : 

1.  construction  by  algebraic  interpolation 

2.  numerical  solution  of  partial  differential  equations 

In  the  latter,  partial  differential  equation  systems  of  elliptic,  parabolic  or 
hyperbolic  type  may  be  used.  Elliptic  systems  may  be  used  to  generate  both 
conformal  and  quasiconformal  mappings,  the  former  being  orthogonal.  In 
this  study,  emphasis  is  on  differential  grid-generation  schemes.  However, 
for  completeness  we  first  describe  algebraic  methods. 

2.1  ALGEBRAIC  METHODS 

The  primary  advantages  of  algebraic  grid  generation  methods  are  that  they 
allow  explicit  control  of  the  physical  grid  shape  and  spacing  and  are  inex¬ 
pensive.  The  fundamental  idea  on  which  most  algebraic  methods  are  based 
is  the  use  of  mathematical  interpolation  functions  to  interpolate  between 
some  known  (or  pre-assigned)  grid  points  (usually  on  the  boundaries),  in  or¬ 
der  to  generate  the  grid  in  between  these  points.  The  interpolation  method 
or  technique  may  vary  with  different  algebraic  methods,  but  the  essential 
idea  remains  the  same.  The  points  between  which  the  interpolation  is  car¬ 
ried  out  need  not  be  boundary  points;  they  can  be  any  specified  points  in 
the  interior  of  the  grid,  through  which  we  desire  that  the  grid  lines  pass. 
The  interpolation  functions  contain  coefficients  which  are  determined  so 
that  the  functional  values  match  the  coordinate  values  at  these  specified 
points.  The  most  difficult  aspect  of  algebraic  grid  generation  is  the  deter¬ 
mination  of  functions  which  control  a  grid. 


For  unidirectional  interpolation  (interpolation  in  one  coordinate  direc¬ 
tion  alone)  the  interpolation  methods  commonly  used  are  polynomial  inter¬ 
polation  and  spline  fitting.  In  polynomial  interpolation  a  single  polynomial 
function  is  used,  to  match  the  coordinate  values  at  all  the  specified  points 
(e.g.,  Langrange  interpolation).  In  some  polynomial  interpolation  meth¬ 
ods  it  is  also  possible  to  specify  the  slope  of  the  function  at  any  or  all  of 
the  specified  grid  points  (e.g.,  Hermite  interpolation).  The  disadvantage 
of  polynomial  interpolation  is  that,  as  the  number  of  specified  quantities 
(i.e.,  points  and/or  slope  at  these  points)  increase,  the  order  of  the  poly¬ 
nomial  increases.  This  may  lead  to  undesirable  oscillations  in  the  behavior 
of  the  polynomial  in  between  the  specified  grid  points.  This  problem  can 
be  avoided  if  spline  fitting  is  used  instead  of  polynomials.  In  spline  fitting, 
instead  of  using  a  single  polynomial  for  all  the  points,  low  order  piecewise 
polynomials  are  used  to  interpolate  the  specified  grid  points.  Slope  con¬ 
tinuity  is  enforced  at  each  point  so  that  smooth  grid  lines  are  obtained. 
The  advantage  of  using  splines  is  that  grid  lines  can  be  generated  relatively 
easily  and  oscillation  of  the  grid  lines  is  reduced.  An  effective  procedure 
based  on  the  description  of  two  exterior  boundaries  and  the  application  of 
either  linear  or  Hermite  cubic  polynomial  interpolation  to  compute  the  in¬ 
terior  grid  is  given  in  Smith  [3].  For  cubic  interpolation  surface  derivatives 
combined  with  magnitude  coefficients  control  the  orthogonality  of  the  grid 
at  and  near  the  boundaries. 

For  example,  non-polynomial  interpolation  functions  are  sometimes  used 
when  the  variation  in  spacing  between  the  specified  grid  points  is  large. 
The  functions  commonly  used  for  this  are  the  exponential  functions  and 
the  hyperbolic  sine  and  tangent  functions.  In  the  multisurface  method,  de¬ 
veloped  by  Eiseman  [4],  interpolation  is  carried  out  between  an  inner  and 
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outer  boundary,  with  a  series  of  control  surfaces  in  between.  This  involves 
the  use  of  field  vectors  tangent  to  the  coordinate  curve  passing  across  these 
surfaces.  The  most  common  form  is  transfinite  interpolation,  which  is  es¬ 
sentially  an  interpolation  between  curves  or  surfaces  rather  than  points. 
It  is  called  “transfinite”  because  it  matches  coordinate  values  on  an  entire 
curve  or  surface.  For  example,  Erikson  [5]  applies  transfinite  interpolation 
for  three  -dimensional  grid  generation  about  wing  -  body  configurations. 
Some  examples  for  airfoils  and  cascade  configurations  are  given  by  Gordon 
and  Thiel  [6]  and  Eiseman  [4,  7]. 

2.2  CONFORMAL  MAPPING  METHODS 

In  conformal  mapping  methods  analytic  functions  are  used  for  transforming 
the  physical  domain  into  intermediate  domains  wherein  the  grid  generation 
problem  is  simpler,  and  subsequently  for  remapping  the  grid  into  the  phys¬ 
ical  or  computational  domain.  These  methods  usually  lead  to  reasonable 
forms  for  the  transformed  partial  differential  equations.  In  addition  they 
can  be  used  for  generating  orthogonal  grids.  A  disadvantage  of  conformal 
methods  is  that  they  provide  little  control  over  the  grid  point  distribution 
for  general  domains. 

The  general  procedure  for  grid  generation  using  conformal  methods  con¬ 
sists  of  two  steps.  The  first  step  consists  of  determining  an  appropriate 
mapping,  or  sequence  of  mappings,  that  would  transform  a  given  physical 
domain  into  a  simple  region.  The  second  step  consists  of  generating  the 
orthogonal  grid  in  the  computational  domain.  This  grid  is  then  remapped 
to  obtain  an  orthogonal  grid  in  the  physical  domain. 

Contours  such  as  single  airfoils  are  generally  mapped  to  “near-circles”  by 
one  or  more  simple  transformations,  and  then  the  “near-circle”  is  mapped 
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to  a  circle  by  a  series  transformation,  (i.e.,  the  Theodorsen  procedure  (Akai 
and  Mueller  [8])).  It  is  necessary  for  convergence  that  “the  near-circle”  be 
sufficiently  close  to  circular.  Also,  there  are  some  very  interesting  ideas 
for  grid  generation  about  multi-component  airfoils  -  Halsey  [9,10]  discusses 
three  methods.  The  first  method  treats  multicomponent  airfoils  by  a  se¬ 
quence  of  transformations  which  map  each  airfoil  to  a  circle  in  succession, 
while  maintaining  previously  established  circles.  The  second  approach  in¬ 
volves  mapping  each  airfoil  to  a  circle  with  no  special  consideration  of 
the  others.  This  process  generally  requires  only  a  few  iterations  to  con¬ 
verge.  The  third  approach  involves  connecting  all  of  the  bodies  in  a  string 
and  mapping  the  resulting  (effective)  single  body.  This  procedure  is  the 
simplest,  but  will  not  give  satisfactory  grids  for  closely  spaced  bodies  in 
general.  Similar  ideas  are  presented  by  Harrington  [11]  where  the  bod¬ 
ies  are  all  mapped  to  rectangles  and  by  Ives  [12]  where  a  Karman-Trefftz 
transformation  is  used  to  map  airfoils  to  the  near  circles. 

Although  the  complex  variable  techniques  by  which  conformal  trans¬ 
formations  are  usually  generated  are  inherently  two-dimensional,  certain 
more  general  cases  can  be  treated  by  rotating  or  stacking  two-dimensional 
systems.  For  example,  Dulikravich  [13]  has  developed  multilevel  three- 
dimensional,  C-type  periodic,  boundary  conforming  grids  for  the  calculation 
of  realistic  turbomachinery  and  propeller  flow  fields.  The  method  is  based 
on  two  analytic  functions  that  conformally  map  a  cascade  of  semi-infinite 
slits  to  a  cascade  of  doubly  infinite  strips  on  different  Riemann  sheets.  For 
the  same  type  of  problems  Anderson,  Davis,  Hankins  and  Edwards  [14]  ap¬ 
ply  Shwarz-  Christoffel  transformation  and  demonstrate  the  validity  of  the 
grid  in  calculations  of  viscous  turbulent  flow  through  a  turbine  inlet  duct. 
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2.3  DIFFERENTIAL  GRID- GENERATION 
SCHEMES 

Differential  grid-generation  schemes  have  received  widespread  attention  be¬ 
cause  of  their  versatility  and  the  ease  with  which  they  can  be  applied.  Here 
the  essential  idea  is  to  generate  the  coordinate  mapping  (and  hence  the 
grid)  as  the  discrete  approximate  solution  to  a  partial  differential  equa¬ 
tion,  matching  the  given  boundary  shape  as  appropriate  boundary  or  initial 
data.  The  underlying  characteristics  of  differential  grid  generation  schemes 
are  depicted  in  fig.  3.  We  first  recall  some  fundamentals  for  second  order 
partial  differential  equations  and  their  solutions. 

2.3.1  GENERAL  PROPERTIES  OF  SECOND  ORDER  PAR¬ 
TIAL  DIFFERENTIAL  EQUATIONS 

The  classification  scheme  for  second  order  PDE’s  depends  on  the  nature  of 
their  characteristics.  In  the  case  of  two  independent  variables,  characteris¬ 
tics  are  lines  in  the  plane  of  the  independent  variables  along  which  “signals” 
can  propagate. 

Consider  the  second  order  quasilinear  PDE  in  two  independent  variables 


of  the  form 


XX  +  B$ xv  ■+■  =  0 


where  A,B,C  and  D  may  all  be  functions  of  x,  y,  $,  and  $y.  The 
classification  of  this  equation  depends  on  the  sign  of  B 2  —  4 AC.  If  B2  — 
4 AC  >  0  the  equation  is  called  hyperbolic;  if  B2  —  4 AC  =  0,  the  equation 
is  parabolic;  if  B2  —  4 AC  <  0,  the  equation  is  elliptic. 

•  HYPERBOLIC  PDE’S 
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Hyperbolic  PDE’s  posses  two  families  of  real  characteristics.  Physical  sys¬ 
tems  that  are  governed  by  hyperbolic  equations  involve  signals  that  prop¬ 
agate  at  finite  speed.  They  are  frequently  posed  in  domains  that  extend 
to  infinity  in  the  timelike  coordinate  and  are  thus  unbounded  in  this  di¬ 
rection.  The  spatial  coordinate  may  or  may  not  be  bounded.  In  either 
case  one  typically  specifies  two  initial  conditions  at  t  =  0.  The  character¬ 
istics  define  the  range  of  influence  of  initial  data.  If  the  spatial  region  is 
bounded,  boundary  conditions  are  also  specified;  otherwise,  we  have  a  pure 
initial-value  problem.  The  prototype  (model)  hyperbolic  equation  is  the 
wave  equation 

-  c2$xx  =  0  (2) 

where  c  is  the  wave  speed. 

•  PARABOLIC  PDE’S 

Parabolic  partial  differential  equations  can  be  regarded  as  the  limit  of  hy¬ 
perbolic  equations  in  which  the  propagation  speed  of  the  signal  becomes 
infinite.  The  prototype  (model)  parabolic  equation  is  the  heat  equation 

=  D*xx  (3) 

Since  equation  (3)  is  first  order  in  time,  only  one  initial  condition  at  t  =  0 
is  necessary,  and  there  is  only  one  set  of  characteristics.  The  domain  of 
solution  is  unbounded  in  time  and  the  spatial  domain  may  be  unbounded 
or  finite. 

•  ELLIPTIC  PDE’S 


Elliptic  partial  differential  equations  have  no  real  characteristics.  In 
elliptic  problems  every  point  in  the  solution  domain  is  effected  by  distur¬ 
bances  at  every  other  point.  The  prototype  (model)  elliptic  PDE  is  the 


Laplace’s  equation 


*xx  +  $  „„  =  0 


Boundary  conditions  are  provided  by  giving  the  value  of  the  dependent 
variable,  its  normal  derivative,  or  some  linear  combination  of  the  two,  at 
each  point  on  the  boundary. 

2.3.2  ITERATIVE  METHODS  FOR  SOLUTION  OF  PDE’S 

Iterative  solutions  and  time  stepping  techniques  to  steady  state  are  closely 
related.  We  consider  iterative  grid  generation  for  elliptic  systems  and 
marching  algorithms  for  parabolic  or  hyperbolic  systems. 

An  iterative  method  applied  to  model  PDE  (4)  is  a  procedure  of  the 
type 

[Ai]$fc+1  =  [A2]$k  +  d  (5) 

where  Ai  and  A2  are  matrices,  d  is  a  vector  and  k  indicates  the  iteration 
level.  Given  some  initial  guess  for  the  solution  ,  we  use  (5)  to  find 
and  <ld1 2 3)  and  so  on.  If  this  is  to  be  a  satisfactory  means  of  solving  the 
original  equation,  it  should  have  the  following  properties: 

1.  it  should  converge  to  the  exact  solution  of  the  original  equation.  That 
is  we  should  have  limfc_00  =  $,  where  $  is  a  solution  of  [A]$  =  b 

2.  the  convergence  should  be  rapid  if  the  method  is  to  be  efficient;  that 
is  the  above  limit  should  be  reached  quickly. 

3.  so  that  each  step  does  not  require  much  computation,  the  matrix  A\ 
should  be  easy  to  invert  and  the  matrix  A2  should  be  as  simple  as 
possible  to  facilitate  the  computation  of  [A]2$*. 


H 


Slvtt*; 


ms 


m 

5 M 


4i.V,y. 

$$ 

m 

If 

*  v!v' 
»»;•*** 

m 


emy: 

tej 


11 


>  t-t  ...I  m  * 


-  IV  fc  1  i 


it 


la  this  section  we  consider  Jacobi,  Gauss-Seidel,  SOR  and  ADI  iterative 
methods 

•  JACOBI  METHOD 

In  the  Jacobi  method  each  new  value  of  a  function  is  computed  entirely 
from  old  values.  Applied  to  model  PDE  equation  (4)  with  uniform  spacing 
Ax  =  A y  =  1  the  method  can  be  written  in  component  form  as 

*{p  =  0.25(S‘+W  +  +  *tn  +  -  0.256„  (6) 

where  i  -  represents  row,  j  -  column  and  k  -  iteration  level.  New  values 
of  are  computed  by  averaging  the  values  of  neighbors  at  the  preceding 
iteration  and  adding  the  inhomogeneous  (boundary  condition)  term. 

•  GAUSS-SEIDEL  METHOD 

When  we  compute  in  a  Jacobi  method,  we  have  already  calculated  the 
values  of  and  $l+ij  .  Since  we  expect  at  each  iteration  level  the  new 

values  to  be  better  approximated  than  the  old  ones,  improvement  could  be 
achieved  by  using  updated  values.  This  leads,  for  model  PDE  eq.  (4),  to 

<!>?+'  =  0.25(*Wj  +  4>f+lj  +  $‘+l,  +  *£„.,)  -  0.256ij  (T) 

which  is  the  basis  of  the  Gauss-Seidel  method.  Sufficient  condition  for 
convergence  of  the  Gauss-Seidel  method  is  diagonal-dominance  of  [A]. 

•  SUCESSIVE  OVER-RELAXATION  (SOR) 

Instead  of  using  1  as  the  input  for  the  next  iteration,  we  can  extrapolate 
the  preceding  results.  This  leads  to  the  procedure  known  as  extrapolation 
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or  over-relaxation.  To  introduce  this  idea,  suppose  that  we  use  Gauss-Seidel 
method  to  compute 

+  dij  (8) 

Using  the  key  idea  of  extrapolation  instead  of  using  as  the  new  iterate 
we  modify  it  by  weighted  averaging  as 


*?!*  =  «*&'  + 


here,  4ft1 2  is  the  most  recent  value  of  $ij  calculated  from  (8),  is  the 
value  from  previous  iteration  and  1  is  the  newly  adjusted  value  of  $,j. 
For  eq.  (9)  to  be  an  extrapolation  and  stable,  overrelaxation  factor  w  should 
be  such  that  1  <  u  <  2.  In  some  problems  (occasionally  for  nonlinear 
problems)  underrelaxation  0  <  u  <  1  is  employed. 


ADI  METHOD 


The  basic  idea  of  ADI  method  is  to  introduce  intermediate  iteration  level 


k  +  1/2  with  operator  splitting.  The  method  can  be  presented  as  follows: 
Let  us  split  matrix  [A]  into  two  matrices  [A]  =  [#]  -f  [V]  such  that 


1.  [H]  +  r[J]  and  [V]  +  r[i]  are  nonsingular  for  any  r  >  0. 


2.  It  is  convenient  to  solve 


([#]  +  r[J])t*i  =  h 


(10a) 


([U]  -I-  r[/])it2  =  b2 


(10b) 


where  bi  and  62  are  any  vectors  and  r  >  0.  Here  the  term  “convenient” 
is  taken  to  mean  ’easy  to  solve  on  a  computer’  and  thus  we  shall  select 
[#]  and  [U]  as  tridiagonal  matrices  (or  matrices  convertible  to  tridiagonal 
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form).  In  the  Peacman-Rachford  method  (Lapidus  and  Pinder  [15]),  we 
select  a  value  r  >  0  and  write 
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([if]  +  r[/])ufc+1/2  =  i  -  ([VI  -  r[I])u‘  (11) 

([V]  +  r[/])u‘+1  =  b  -  ([H]  -  r[/])u*+V!  (12) 

with  the  vector  nfc+1/2  being  an  intermediate  vector.  Thus  (11)  is  used  with 
the  matrices  [H]  and  [V],  and  r  to  calculate  u^+1^2);  then  (12)  is  used 
with  r  to  evaluate  u^+1^.  In  this  way  the  iteration  proceeds  from  u^°\  u M... 
.  The  idea  here  is  the  selection  of  [H]  and  [V]  so  that  (11)  uses  lines  in 
the  x  direction  only  and  (12)  uses  lines  in  the  y  direction  only.  In  order  to 
consider  convergence  of  the  method  let  us  define  error  vector  —  u, 

where  u,  the  correct  solution,  satisfies 


It  then  follows  that 


m + mu = b 


e(fc+D  =  [£],«<*> 


where  [G]r  =  ([F]4-J'[/])_1([jfT]— r[/])([.£f]+r[/])-1([K]  — r[I])  is  the  iteration 
matrix.  It  is  necessary  for  convergence  that  the  spectral  radius  [G]r  satisfies 
r([G]r)  <  1.  If  matrices  [H]  and  [V]  are  positive  define,  the  above  condition 
holds.  For  model  equation  (4)  the  method  can  be  written  (with  r  =  1)  as 


first  step: 


4Ci'S  -  Ml?'2  +  <C'.?  =  -Mlu  -  M‘j  +  +Uu)  04) 


second  step: 


-  Mt?  +  Mtj  =  -Mi??  -  Mi?'2 + *is?) 
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The  iteration  begin  with  an  initial  estimate  of  the  $1°)  vector,  and  com¬ 
putation  is  terminated  when  successive  calculations  agree  within  a  given 
tolerance. 

Except  for  the  additional  work  of  recalculating  right  sides,  this  method 
has  little  computational  effort  and  converges  fast.  The  coefficients  matrices 
are  always  the  same,  so  the  reduction  step  need  be  performed  only  once. 


2.4  ELLIPTIC  GRID  GENERATION  SCHEMES 

Elliptic  schemes  require  the  specification  of  data  on  the  entire  boundary  of 
the  domain.  The  location  of  the  interior  grid  points  is  then  determined  by 
solving  a  set  of  elliptic  partial  differential  equations.  The  most  commonly 
used  are  the  Laplace  and  Poisson  equations.  In  the  former,  with  equispread 
boundary  points,  the  coordinate  lines  will  tend  to  be  equally  spaced  in 
the  absence  of  boundary  curvature,  because  of  the  smoothing  effect  of  the 
Laplacian.  In  the  Poisson  case  mesh  control  (coordinate  lines  clustering 
and  inclination)  can  be  accomplished  by  means  of  an  appropriate  choice  of 
forcing  functions. 

The  Poisson  system  defines  the  computational  (£,  77)  coordinates  by 


£xx  4"  £yy  —  1} ) 

Tfcx  4  *7yy  —  7) 


(16a) 

(16b) 


where  P(£ ,  77)  and  Q(£,  77)  represent  forcing  functions  that  are  specified  and 
depend  on  the  nature  of  the  problem  and  the  desired  grid  point  distribution. 

The  primary  advantages  of  such  elliptic  systems  for  grid  generation 
are:  (1)  they  frequently  generate  one-to-one  mapping;  (2)  the  coordinate 
lines  are  smooth  owtng  to  the  inherent  smoothness  in  the  solution  of  elliptic 
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systems;  and  (3)  boundary  slope  discontinuities  are  not  propagated  into  the 
field. 

Grid  generation  schemes  formulated  using  the  Laplace  and  Poisson 
equations  have  been  extensively  applied  for  various  two-  and  three-dimensio¬ 
nal  configurations.  For  example,  Thompson,  Thames  and  Mastin  [16]  and 
Holst  [17]  control  the  interior  grid  distribution  using  exponential  forcing 
functions  which  contain  adjustable  parameters.  The  selection  of  these  pa¬ 
rameters,  however,  depends  on  the  shape  of  the  body  surface.  They  provide 
the  means  to  move  the  coordinate  lines  around  but  not  the  proper  ampli¬ 
tudes  and  decay  factors  necessary  to  achieve  desired  spacing  distributions. 
Hence,  this  technique  requires  interaction  from  the  user. 

Sorenson  and  Steger  [18]  introduce  automatic  mesh-point  clustering 
near  the  boundary  and  Sorenson  [19],  applies  the  idea  to  elliptic  grid  gen¬ 
eration  for  discretizing  the  flow  about  the  augmentor-  wing  using  a  compo¬ 
nent  adaptive  grid  interfacing  technique.  His  method  differs  from  that  of 
Thompson,  Thames,  and  Mastin  [16]  in  that  it  uses  forcing  terms  P{£,tj) 
and  Q{£,  77),  which  yield  the  ability  to  arbitrarily  impose  two  types  of  con¬ 
trol  on  the  grid  at  boundaries.  The  first  type  of  control  is  of  the  angle 
on  inclination  at  which  the  lines  of  constant  £  intersect  the  boundaries; 
the  second  type  of  control  is  on  the  distance  between  the  airfoil  and  lines 
=  constant. 

Coleman  [20]  applies  the  Laplace  or  Poisson  equations  to  produce  a 
mesh  which  adapts  to  the  physical  region  of  interest.  The  method  provides 
arbitrary  segmentation  in  the  physical  plane  and  mapping  to  a  union  of 
rectangular  grid  sections. 

Chen  and  Obasih  [21]  introduce  contracting  functions  for  the  grid  inside 
the  computational  domain.  They  treat  contraction  functions  P(£,tj)  and 
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I  Q(£,r))  iQ  the  original  Poisson  equation  as  unknowns,  and  require  orthog- 

|  onality  of  coordinate  lines  near  boundaries,  to  determine  expressions  for 

'  P(£,t;)  an<t  Q(£,77)-  Similar  ideas  are  given  by  Visbal  and  Knight  [22]  and 

j  Shieh  [23]. 

2.5  HYPERBOLIC  SCHEMES 

Hyperbolic  schemes  are  obtained  by  integrating  a  hyperbolic  PDE  for  the 
grid  lines  and  require  Cauchy  data  on  the  boundary.  The  grid  may  be 
generated  explicitly  by  marching  outward  from  the  initial  boundary.  For 
example  Thompson,  Warsi  and  Mastin  [16]  discuss  the  system 

^Vv-y^v  =  v(z,t))  (n) 

xixv  -  VtVv  =  0  (18) 

for  generating  a  grid.  Here  V(£,  rj)  is  the  Jacobian  of  the  transformation 
and  represents  the  area  in  physical  space  for  a  given  unit  area  element 
in  computational  space.  If  ^(£,77)  is  given  as  function  of  position  then 
equation  (17)  can  be  used  for  grid  spacing  control  in  the  physical  plane. 
Equation  (18)  represents  a  measure  of  the  orthogonality  of  grid  lines  in  the 
physical  space.  After  linearization  of  equations  (17,18)  about  a  known  state 
( x,y )  we  have  the  system 

[*A]r*  +  [B]rv  =  /  (19) 

where 


X 

.  y . 

;M  = 

Vf, 

.  Vr]  ~xv  . 

;[s]  = 

H 

Vv  ' 

;/  = 

0 

_  K+  V  _ 

The  eigenvalues  of  [B\  1  [A] 


are  real  and  positive  which  implies  that  system  (19)  is  hyperbolic  in  the  77 
direction  and  can  be  marched  in  77  as  long  as  x%  +  y\  ^  0 

In  generating  a  grid  with  this  scheme  first  assume  the  body  surface  is 
the  77  =  0  surface  and  specify  the  distribution  of  points  along  the  body. 
Next  the  quantity  V(£,r])  in  eq.(17)  is  required.  Steger  and  Sorenson  [25] 
suggest  that  V(£,t))  can  be  determined  by  laying  out  a  straight  line  with 
length  equal  to  that  of  the  body  surface,  and  then  lay  out  the  body  point 
distribution  on  this  line  (i.e.,  arc  length  distribution  of  surface  points). 
Next  a  line  parallel  to  the  first  line  is  drawn  on  the  77  =  constant  surface  as 
desired.  Once  this  is  done,  the  quantity  V(£,  77)  is  determined  by  estimating 
the  area  elements  of  the  grid.  Application  of  a  hyperbolic  scheme  for  airfoil 
problems  is  given  by  Steger  and  Chaussee  [26]  and  for  the  three-dimensional 
grid  about  a  shuttle  cross-section  by  Kutler  [27]. 
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2.6  PARABOLIC  SCHEME 

Grids  may  be  generated  by  parabolic  schemes  explicitly  by  marching  as 
in  the  hyperbolic  grid-generation  method  from  an  “initial”  boundary.  As 
in  the  elliptic  schemes  diffusion  smoothes  out  any  irregularities  in  the  grid 
due  to  the  shape  of  the  initial  line  (inner  boundary).  The  importance  of  the 
marching  algorithm  is  twofold:  first,  the  computational  time  required  may 
frequently  be  only  a  very  small  fraction  of  that  for  elliptic  grid  generation ; 
second,  storage  required  during  grid  generation  can  be  substantially  reduced 
from  that  required  by  the  elliptic  grid  generation  method.  We  demonstrate 
these  points  for  a  particular  form  of  parabolic  scheme  in  later  numerical 
experiments. 

The  parabolic  schemes  are  applicable  to  both,  two-  and  three-dimensi¬ 
onal  problems.  Nakamura  [28]  applies  a  parabolic  scheme  with  simplified 
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second  derivatives  to  O-type  and  H-type  grids  for  single  airfoils  and  Naka¬ 
mura  [29]  applies  parabolic  differencing  scheme  to  one  coordinate  variable 
of  a  three-dimensional  elliptic  equation  set  for  fuselage-wing  configurations. 
Edwards  [30]  extends  this  three-dimensional  grid  generation  algorithm  to 
two  of  the  coordinate  variables  while  the  third  is  centrally  differenced  and 
generates  a  grid  about  a  wing-fuselage  and  an  aircraft  configuration. 

Parabolic  systems  have  previously  been  considered  to  be  limited  in  the 
configurations  that  they  can  handle.  In  the  present  study  we  demonstrate 
that  parabolic  schemes  can  be  extended  to  complex  configurations  with 
orthogonality  satisfied  on  inner  boundaries,  and  apply  the  method  to  rep¬ 
resentative  multi-airfoil  problems. 

2.7  ADAPTIVE  GRIDS 

Adaptive  grids  are  dynamic  grids  in  which  the  grid  points  are  automatically 
readjusted  as  the  solution  evolves.  Some  aspect  of  the  developing  solution 
must  be  used  as  the  error  indicator  for  redistributing  the  grid  points.  We 
note  that  readjustment  of  the  grid  points  must  yield  a  grid  of  acceptable 
smoothness  and  orthogonality.  The  grid  points  should  be  clustered  in  re¬ 
gions  of  large  solution  variations  and  be  well  graded  into  such  regions.  As 
an  example,  Brackbill  [31]  employs  an  objective  function  which  contains 
a  measure  of  grid  smoothness,  orthogonality,  and  volume  variation  that  is 
minimized  over  a  class  of  admissible  grids.  The  smoothness  of  the  trans¬ 
formation  is  represented  by  the  integral 


Is  =  /J(V£)2  +  (Vr,f\iV 

JD 


A  measure  of  orthogonality  is  provided  by 


Io=  f  (V^rj)2J3dV 
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and  an  error  indicator  measure  (weighted  by  the  Jacobian) 


Iv  ~  / 
Jd 


WJdV 


where  W  is  a  given  weighting  function  and  J  is  a  measure  of  cell  vol¬ 
ume,  usually  the  Jacobian  of  transformation.  If  Iy  is  minimized  for  a  fixed 
W(;c,y),  J  (the  grid  size)  is  made  small  where  W  is  large  and  vice  versa. 
The  weight  function  W(x,y)  is  the  error  indicator  and  is  to  be  a  function 
of  some  measure  of  the  solution  error,  so  that  the  spacing  will  be  reduced 
where  the  error  is  large.  As  a  final  objective  function  a  weighted  sum  of 
(21)  -  (23)  is  taken.  The  inclusion  of  (21)  and  (22)  as  integral  constraints  in 
the  objective  function  ensures  that  the  optimal  grid  has  reasonable  smooth¬ 
ness  and  orthogonality,  respectively.  Brackbill  [31]  implements  the  adap¬ 
tive  control  through  terms  in  an  elliptic  generating  system  and  Saltzman 
and  Brackbill  [32]  demonstrate  results  for  multiple  shock  reflections  in  a 
wind  tunnel  problem.  Anderson  and  Rai  [33]  give  a  procedure  based  on 
an  analogy  with  electrostatic  charge  attraction  which  is  applicable  to  any 
coordinate  system  and  demonstrate  the  effectiveness  of  the  shock  aligning 
scheme  for  a  straight  oblique  shock  in  a  uniform  supersonic  flow. 

3  NEED  FOR  BOUNDARY-FITTED  CO¬ 
ORDINATE  SYSTEM 

There  arises  in  all  problems  concerned  with  the  numerical  solution  of  par¬ 
tial  differential  equations  the  need  for  accurate  numerical  representation  of 
boundary  conditions.  For  example,  in  the  present  study  we  seek  to  define 
points  of  a  finite  difference  grid  constructed  on  coordinate  lines  that  coin¬ 
cide  with  the  boundary.  That  is  one  coordinate  variable  can  be  specified  to 
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be  constant  on  each  of  boundaries,  and  a  monotonic  variation  of  the  other 
coordinate  around  each  boundary  can  be  specified  as  the  data  for  our  dis¬ 
cretized  partial  differential  equation  problem.  It  then  remains  to  generate 
values  of  these  coordinates  in  the  field  from  prescribed  boundary  or  initial 
values.  There  must,  of  course,  be  a  unique  correspondence  between  the 
basic  coordinate  system  and  the  curvilinear  coordinates;  i.e.,  the  mapping 
of  the  physical  region  onto  a  transformed  computational  region  must  be 
one-to-one,  so  that  every  point  in  the  physical  field  corresponds  to  one, 
and  only  one,  point  in  the  transformed  field,  and  vice  versa.  Coordinate 
lines  of  the  same  family  must  not  cross,  and  lines  of  different  families  must 
not  cross  more  than  once  (i.e.,  the  Jacobian  of  the  transformation  must  be 
nonsingular  at  each  point).  Further,  the  grid  lines  should  be  smooth  to  pro¬ 
vide  continuous  transformation  derivatives.  Grid  points  should  be  closely 
spaced  in  the  physical  domain  where  large  numerical  errors  are  expected 
and  excessive  grid  skewness  should  be  avoided  because  it  can  sometimes 
amplify  truncation  error. 

Since  the  boundary-fitted  coordinate  system  has  coordinate  lines  coin¬ 
cident  with  the  surface  contours  of  all  bodies  present,  all  boundary  condi¬ 
tions  for  the  problem  can  be  expressed  at  grid  points  without  recourse  to 
interpolation  or  projection,  and  normal  derivatives  on  the  bodies  can  be 
represented  using  finite  differences  between  grid  points  even  though  the  co¬ 
ordinate  system  may  not  be  orthogonal  at  the  boundary.  The  transformed 
equations  can  then  be  approximated  using  finite  difference  expressions  and 
solved  numerically  in  the  transformed  plane. 
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3.1  TYPE  OF  GRIDS 

Requirements  on  the  type  of  grid  to  be  used  in  the  physical  domain  are 
mainly  problem  dependent.  Three  grid  types  for  aerodynamic  flow  prob¬ 
lems  are  in  common  use  and  have  been  termed  H-type,  O-type  and  C-type. 

•  A  H-type  grid  (fig.  4a)  provides  excellent  resolution  of  the  flow  field 
at  upstream  and  downstream  infinity.  It  is  also  the  simplest  grid  to 
generate.  At  the  same  time,  H-type  grids  do  not  provide  an  accurate 
treatment  of  rounded  leading  and  trailing  edges  and  grid  points  in 
the  flow  domain  away  from  the  boundaries  are  not  well  clustered. 

•  An  O-type  grid  (fig.  4b)  represents  a  coordinate  system  having  lines 
encircling  a  body.  It  gives  very  poor  resolution  at  infinity,  but  pro¬ 
vides  very  good  resolution  for  blunt  and  rounded  edges,  and  uses  few 
grid  points 

•  A  C-type  grid  (fig.  4c)  indicates  a  coordinate  system  with  lines  ema¬ 
nating  from  a  boundary,  passing  around  a  body  and  returning  to  the 
boundary.  The  grid  represents  a  combination  of  an  O-type  grid  in 
the  upstream  region  and  H-type  grid  in  the  downstream  region.  This 
type  of  grid  provides  a  good  treatment  of  all  boundary  and  period¬ 
icity  conditions  including  wake  treatment  and  supersonic  exit  flow, 
although  it  may  not  yield  an  adequate  resolution  at  upstream  infinity 
for  certain  applications  such  as  cascades. 
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3.2  TRANSFORMATIONS  BETWEEN  PHYSICAL 
AND  COMPUTATIONAL  PLANE 

The  general  transformation  from  the  physical  plane  [ x ,  j/]  to  the  transformed 
plane  has  the  form 

0  =  U(^)1  (24\ 

t?  rj(x,y)  ' 


The  Jacobian  matrix  for  this  transformation  is 


-  <MiHl  ~  @  1  _ 

d(x,y)  1  fe  it  J 


r]x  r)y 


The  inverse  function  or  transformation  of  (24)  is, 


x  1  _  x[Z,ri) 

y  [y(Z,n) ' 


with  Jacobian  matrix 


_  sm  _  r  f  f  i  _  r  *<  *. 


The  Jacobian  is  then: 


Since  J\  =  J2-1, 


^(C.7?)  L  If  J  L  &  Vo 


J  =  det[J ]  =  -  Xr,y( 


£x  £y  __  j- 1  2/>7  ~~XV 

Vr  Vv  .  % 


we  get  the  following  relations: 
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(30b) 


Using  the  chain  rule,  partial  derivatives  of  a  given  function  /  with  re¬ 
spect  to  x  and  y  are  transformed  as  follows, 

f  _  Of  _  Jj&S  _  Vnh  -  Vjfn  \ 

Jx  ~  dx  ~  ”  j  w 

Of  &(xif)  _  r  _  f 

f  _  ££  _  /«in\ 

Iv~  dy~  Steal  “  J  1  j 

fxx  dx(  j  ’ 

=  (i/J/«  -  Zvivnh*  +  y?/w)/ /2 

+  [(y?y«  -  2«ejfc,y{,  +  y?yw)(z„/«  -  *«/») 

+  [(»?*«  -  2 y*y„y*„  +  ylxm)(ydv  -  yJt)\/J3  (33) 

/w  ~  J  ’ 

=  (*J/«  "  2xi^hr,  +  X\fm)/J 2 

"b  [(*,y«  —  ^xixvyiv  "b  xtyvv)(xvfi  ~  xzfv) 

+  l(*J*«  -  2a**„*?„  +  z|*™)(y«/»  -  yn/e)]/*^3  (34) 

/xy  =  \(.x(yri  "b  *vVt)ftn  —  x^y^fw )/J 

+  [*f»y»?*«  -  (Xtyr,  +  *nVl)*ln  +  xmxrm\{yr,h  -  ytfv)/J3 
+  lxvyr,yu  -  (*<y»  +  *»y«)y^  +  xmvm )Mv  -  xvh)N3  (35) 
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and  so  on. 

Sufficient  conditions  for  the  transformations  (26)  to  exist  are  given  by 
the  inverse  function  theorem  that  states:  if  the  component  functions  of  (24) 
are  continuously  differentiable  at  some  point,  say  (a?i,2/i),  and  Jacobian 
matrix  (25)  is  nonsingular  at  ( a?i ,  3/1 )  then  there  exists  a  region  M\  about 
(zi,2/i)  such  that  the  inverse  function  (26)  exists  and  (28)  holds  for  all 
[(z:,j/)]  in  M\.  It  is  apparent  that  the  theorem  guarantees  existence  only 
in  a  local  sense.  For  this  reason  component  functions  of  (24)  which  posses 
even  more  desirable  properties  than  those  stated  for  the  inverse  function 
theorem  are  sought. 

4  AN  ELLIPTIC  SCHEME  AND 
APPLICATION 

The  basic  idea  of  the  transformations  in  section  (3)  is  to  let  the  compo¬ 
nent  functions  of  (24)  be  solutions  of  an  elliptic  Dirichlet  boundary  value 
problem.  An  obvious  choice  is  to  require  that  f(ar,  y)  and  r)(x,y )  be  ei¬ 
ther  harmonic,  subharmonic,  or  superharmonic.  Harmonic  functions  obey 
a  maximum  principle,  which  states  that  the  maximum  and  minimum  val¬ 
ues  of  the  function  must  occur  on  the  boundaries  of  the  region  D.  Since 
no  extremes  occur  within  D,  the  first  derivatives  of  the  function  will  not 
simultaneously  vanishes  in  D,  and  hence  the  Jacobian  will  not  be  zero  due 
to  the  presence  of  an  extremum.  The  maximum  principle  also  guarantees 
uniqueness  of  the  coordinate  functions  £(x,  y)  and  rj(x,y),  and  thus  ensures 
that  no  overlapping  of  the  boundaries  will  occur. 

Now,  let  us  take  Laplace’s  equation  as  the  generating  elliptic  system. 


£xx  "b  £vv  —  ® 
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*]xx  +  TJyy  =  0 

with  the  Dirichlet  boundary  conditions 


6(*,y) 

.  *1  . 
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’  6(ar,y) 
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(36b) 

(37a) 

(37b) 


where  rfr  and  %  are  constants  and  £1(2:,  y)  and  £2(2;,  y )  are  specified  mono¬ 
tonic  functions  on  G\  and  Gi  respectively  (fig.  6).  That  is,  let  £(x,  y) 
and  rj(x,y )  be  harmonic  in  D.  This  generation  system  guarantees  a  one- 
to-one  mapping  for  boundary-conforming  curvilinear  coordinate  systems 
on  general  closed  boundaries.  Since  all  numerical  computations  are  to  be 
performed  in  a  uniform  rectangular  transformed  plane,  the  dependent  and 
independent  variables  must  be  interchanged  in  (36).  Using  transformations 
(31)  -  (35),  and  knowing  that  coordinate  lines  in  the  transformed  plane  are 
constant,  we  get 

£■*  =  i(y$yx  -  ZytyvVtv  +  vhvn)xv 

+  (yWi  -  2 ytvivir,  +  y\xrm){-yr,)V j3  (38) 

Ixx  =  -  ZvtlhVto  +  J^Iftw)(-*{) 

+  (ylxtt  -  tytunv* »  +  ylxw)(yc))/J3  (39) 

£yy  =  —  2x(xvl Kv  *b  x^yvv)xv 

*b  i.xr)Xii  ~  2xZxvxtv  "b  x£xvv)(~ Vv)\l ^  (40) 

Vyy  =  [(*Jy«  ~  2 +  xJjfcw)(-*{) 

"b  (x„x«  —  2x(XvX(r)  ■+■  x^xm)(y{)\/ J  (41) 
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and  the  transformed  system  has  the  following  form 

(Axu  +  Bxir)  +  Cx^i-Vr,)  4-  ( Ayu  +  Byiv  4-  Cyrm)(xr))  =  0  (42a) 

(Axu  +  Bxtr,  +  Cx„„)(y?)  4-  (Ayu  +  Bytr)  +  Cy„„)(-x€)  =  0  (42b) 


Where  the  coefficients  are 


A  =  xl  +  yl 


B  =  —2(x^xv  4-  y^r,) 
C=x\  +  y\ 


—  C\yv  +  C2x„  =  0 

Ciyt  -  C2x*  =  0 


(43a) 

(43b) 

(43c) 


(44a) 

(44b) 


where 


Ci  —  Ax^£  4"  Bx^rj  4"  Cxnn 

C2  =  Ay^i  4-  By^r,  +  Cyvv 

Since  the  Jacobian  determinant  is  nonzero,  a  nontrivial  solution  of  system 
(44)  exists  if  Ci  =  C2  =  0,  and  we  get  the  coupled  system 


Ax^^  4"  Bx^  4"  Cxjpj  —  0 

+  By  &  +  Cyvv  =  o 


with  transformed  boundary  conditions 


*  1  _  I"  MZ,m)  1  . 
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(45b) 
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The  functions  Pi(£,  ?7i),p2(£,  ??i),  tfi(f,  772)  and  ?2(^!%)  are  specified  by  the 
known  shape  of  the  inner  and  outer  boundary.  Assuming  uniform  spacing 
on  the  computational  domain,  the  central  difference  equations  approximat¬ 
ing  (45)  at  grid  point  (i,j)  may  be  written, 

A(ar,_i  j  —  2xitj  4-  £«+i,j)  + 

~  xi- ij+i  —  £«+ lj-i  +  *»+ij+i)/4  + 

C(xi,j- i  —  2x ij  +  r,,j+i)  =  0  (47) 

~  4-  jfc+ij)  4- 

•S(y.-i,j-i  —  —  y,+i,j_i  +  yi+1j+1)/4  + 

C(yi,i-1  -  2yij  +  y,j+i)  =  0  (48) 

that  is,  in  general  form 

—  2«,j  4-  s.+ij)  + 

~  ft-lj+l  ~  6i+lJ-l  4  Si+ij+i)/4  -f 

—  2&ij  4-  Sij+i)  =  0  (49) 

where  s,j  represents  Cartesian  coordinate  x  or  y. 

The  discretized  system  (49)  can  be  solved  iteratively  using  a  relaxation 
scheme  of  the  form  described  previously,  such  as  successive  over- relaxation 
(SOR)  or  alternating  direction  implicit  (ADI),  until  a  specified  convergence 
criterion  is  satisfied. 

Adding  forcing  functions  terms  in  eqs.  (36)  we  get  the  Poisson  system 


4-  —  P{£,  y) 


(50a) 
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Vxx  +  T)yy  =  Q(Z,ri) 


(50b) 


and  applying  the  same  transformation  procedure  this  yields  instead  of 
(42) 

AX(t  +  Bx(„  +  Cx„„  =  -JJ[.P( {,t?)*f  +  Q((,  (51a) 


-4»i  +  By(,  +  Cym  =  -JJ[P({.  r/)y(  +  Q({, 


(51b) 


As  noted  previously,  the  distribution  of  grid  lines  obtained  by  solution 
of  the  discretized  form  of  equations  (49)  will  tend  to  be  equally  spaced  in 
the  absence  of  boundary  curvature  because  of  the  smoothing  effect  of  the 
Laplacian,  but  will  become  more  closely  spaced  over  convex  boundaries, 
and  less  so  over  concave  boundaries,  as  illustrated  in  fig.  5a.  For  the  first 
example  in  fig.  5a  we  have  rjxx  >  0  because  of  the  convex  curvature  of 
the  lines  of  constant  r).  Therefore  it  follows  that  r)yy  <  0,  and  the  spacing 
between  the  77  -  lines  must  increase  with  y.  The  77  -  lines  thus  will  tend  to 
be  more  closely  spaced  over  a  such  convex  boundary  segment.  For  concave 
regions  (the  second  example  in  fig.  5a)  we  have  T]xx  <  0,  so  that  T)yv  must 
be  positive,  and  hence  the  spacing  of  the  77  -  lines  must  decrease  outward 
from  this  concave  boundary.  Considering  eq.  (50b),  it  follows  that  negative 
values  of  Q(£,  77)  will  tend  to  cause  the  coordinate  line  spacing  to  increase 
more  rapidly  outward  from  the  boundary.  In  general,  negative  values  of 
the  control  function  Q(£,  77)  will  cause  the  77  -  lines  to  tend  to  move  in  the 
direction  of  decreasing  77,  while  negative  values  of  P(£,i 7)  in  eq.  (50a)  will 
cause  £  -  lines  to  tend  to  move  in  the  direction  of  decreasing  £.  These  effects 
are  illustrated  in  fig.  5b  for  an  77  -  line  boundary.  With  the  boundary  values 
fixed,  the  £  -  lines  here  cannot  change  the  intersection  with  the  boundary. 
The  effect  of  control  function  P(£,r})  in  this  case  is  to  change  the  angle  of 
intersection  at  the  boundary,  causing  £  -  lines  to  “lean”  in  the  direction  of 
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decreasing  £. 

Generalizing,  a  negative  value  of  the  Laplacian  of  one  of  the  curvilinear 
coordinates  causes  the  lines  on  which  that  coordinate  is  constant  to  move 
in  the  direction  in  which  that  coordinate  decreases.  Positive  values  of  the 
Laplacian  naturally  result  in  the  opposite  effect. 

4.1  COORDINATE  SYSTEM  CONTROL 

Control  of  the  spacing  of  the  coordinate  lines  on  the  body  is  accomplished 
through  the  input  data,  while  the  spacing  of  the  coordinate  lines  in  the  field 
is  controlled  by  varying  the  elliptic  generating  system  for  the  coordinates. 
One  procedure  that  has  proven  effective  for  some  aerodynamic  applications 
(  Holst  [17],  Thompson,  Thames  and  Mastin  [16]  ),  is  to  choose  forcing 
functions  P(£,t])  and  Q(£,rj)  as  exponential  terms 

PU,rj)  =  -  ^2  aisign(£  -  &)e(-c,K~'e,|) 

t=i 

M  - 

—  bjsign(£  —  £j)e~d]  V^-l>)a+(T>-»h)a}  (52) 

j=i 


Q(Z,V)  =  -X]aisisfn(77-7?1)e(  C,,T'  *l) 

«=i 

M  - 

—  S  bjsign(rj  —  77i)e-d;V^-*-d2+(,!-T7j)2] 


The  first  terms  in  P(£,t))  and  Q{£,r})  have  the  effect  of  attracting  the  £  = 
constant  lines  to  the  £  =  £,  lines  in  eq.  (52),  and  attracting  77  =  constant 
lines  to  77  =  77*  lines  in  eq.  (53).  The  second  terms  cause  £  =  constant  lines 
to  be  attracted  to  the  points  (£;  ,  77;)  in  eq.  (52)  with  a  similar  effect  on 
77  =  constant  lines  in  eq  (53).  Introducing  the  sign  function,  the  equations 
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(52)  and  (53)  are  no  longer  subharmonic  or  superharmonic,  since  the  sign 
function  causes  a  sign  change  on  the  right  side  when  the  attraction  is  to 
lines  or  points  not  on  the  boundaries.  So,  it  is  possible  that  too  strong  an 
attraction  may  cause  the  system  to  overlap  and  produce  unusable  grids. 
The  use  of  the  sign  -  changing  function  is  used  to  cause  attraction  to  both 
sides  of  a  line  or  point  in  the  field.  Elimination  of  this  function  causes 
attraction  on  one  side  and  repulsion  on  the  other.  If  it  is  only  desired  that 
coordinate  lines  be  concentrated  near  one  boundary,  (i.e.,  the  body  surface) 
then  there  is  no  need  for  the  sign  change,  and  the  sign  function  can  be 
eliminated.  In  this  case  the  equations  are  subharmonic  or  superharmonic, 
and  a  maximum  principle  is  in  effect  to  prevent  overlap.  The  effect  of  the 
amplitude  factors  in  (52,  53)  is  shown  in  fig.  5c. 

4.2  APPLICATION  OF  THE  ELLIPTIC  SCHEME 

The  elliptic  grid  generation  scheme  can  be  applied  to  both  single  and  mul¬ 
tiple  bodies.  Referring  to  fig.  6,  assume  that  the  body  contour  and  the 
outer  boundary  transform,  respectively,  to  the  rj  -  lines  forming  the  lower 
and  upper  sides  of  the  rectangular  transformed  plane;  arbitrary  cuts  that 
join  those  boundaries  map  to  the  £  -  lines  forming  left  and  right  sides  of 
the  transformed  plane.  Thus  the  left  and  right  side  vertical  boundaries  in 
the  computational  plane  are  coincident  in  the  physical  plane;  the  values  of 
x  and  y  coordinates  are  equal  along  these  lines.  The  computational  field 
size  is  (J M  —  2)  x  ( IM  —  1).  Boundary  values  are  specified  on  j  =  1  and 
j  —  JM  for  all  1  <  i  <  IM.  The  line  j  =  1  corresponds  to  the  body  sur¬ 
face  in  the  physical  plane  while  j  =  JM  corresponds  to  the  outer  boundary. 
The  discretized  form  of  the  governing  equation  in  the  computational  plane 
is  given  by  (49).  The  re  -  entrant  boundaries  occur  at  i  =  1  and  i  =  IM. 
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Since  the  values  of  x  and  y  are  equal  along  these  lines,  iteration  is  necessary 
along  only  one  of  them.  For  i  =  1,  the  £  -  derivatives  along  this  line  can  be 
approximated  as 

=  (*2 j  -  *imj)I2  (54a) 

=  x2,j  —  2xi,j  +  xim- i,j  (54b) 

(x(T>)iJ  =  (X2J+1  —  X2J-1  +  XIM-1J-1  -  xIM-l,j+l)/4  (54c) 


for  2  <  j  <  JM—l.  Similar  expressions  can  be  used  for  the  derivatives  of  y. 
The  set  of  simultaneous  difference  equations  in  x and  is  solved  by  point 
SOR  iteration  in  the  present  numerical  experiments.  The  same  idea  may 
be  extended  to  multiple  bodies.  In  this  case  we  introduce  additional  cuts 
that  connect  bodies  and  apply  a  computational  procedure  similar  to  that 
for  the  single  body.  For  two-component  bodies  the  transformed  boundary 
conditions  are  given  as  follows  (fig.  7) 
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The  functions  p1,p2,  9i,  72,  P3,P4, 93  and  94  are  specified  by  the  known  shape 
of  the  contours  G\,  G2,  G7  and  G8  and  the  specified  distribution  of  £  thereon. 
Re-entrant  boundaries  occur  on  G l,  G\,  G\  and  Gg.  Derivative  approxima¬ 
tions  on  these  boundaries  are  determined  using  a  procedure  similar  to  that 
given  by  (54).  In  the  present  study  we  generate  grids  about  single  and 


multiple  airfoils  by  use  of  equations  (45)  and  (51)  on  uniform  grids  in  the 
computational  plane  following  the  procedure  of  Thompson,  Thames  and 
Mastin  [16]. 

Fig.  8  shows  an  O-grid  about  NACA  0012  airfoil  generated  by  the  ellip¬ 
tic  scheme  with  clustering  functions  P(£,r))  and  Q(£,r?)  given  by  eqs.(52, 
53).  The  grid  size  is  61  x  28  and  radius  of  the  outer  boundary  3.0  x  chord. 
The  attraction  is  applied  for  five  coordinate  lines  near  the  airfoil  surface 
and  leading  and  trailing  edges.  The  initial  guess  is  determined  using  the 
average  values  of  four  boundary  points.  For  a  fixed  acceleration  parameter, 
a  convergent  solution  was  obtained  after  99  iterations  and  48.561  seconds 
of  CPU  time  on  the  CDC  750/175.  The  grid  in  fig.  9  is  generated  using 
the  elliptic  scheme  (45).  The  grid  size  is  49  x  16  and  radius  of  the  outer 
boundary  1.2  x  chord.  The  initial  guess  is  determined  as  the  solution  of  the 
parabolic  scheme  (68,  69).  In  this  case  a  convergent  solution  was  obtained 
in  1 .971  seconds  of  CPU  time.  Good  clustering  of  the  coordinate  lines  at 
leading  and  trailing  edges  of  the  airfoil  is  evident. 

In  fig.  10  and  fig.  11  C-type  and  H-type  grids,  respectively,  are  shown. 
Both  grids  are  generated  using  the  elliptic  system  (45).  The  initial  guess 
is  determined  in  the  same  way  as  for  the  O-type  grid.  It  can  be  seen  in 
all  three  cases  that  the  system  (45)  provides  good  smoothness  of  the  grid 
lines,  but  does  not  provide  clustering  of  the  grid  lines  at  the  body  surface. 

Fig.  12  is  an  O-grid  about  two  NACA  0012  airfoils  with  flap  angle 
of  25  degrees.  Grid  size  is  69  x  20  and  radius  of  the  outer  boundary  4.0 
x  chord.  Attraction  is  applied  for  six  coordinate  lines  near  the  airfoils 
and  at  the  leading  and  trailing  edges  of  both  airfoils.  The  initial  guess  is 
specified  using  the  average  values  of  four  boundary  points.  Convergence 
of  the  solution,  with  constant  acceleration  factor,  was  obtained  after  275 
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iterations  and  105.702  seconds  of  CPU  time.  A  detailed  grid  distribution 
near  the  airfoils  is  given  in  fig.  13. 

In  all  numerical  experiments  it  was  observed  that  the  efficiency  of  the 
elliptic  scheme  depends  on  the  accuracy  of  the  initial  guess.  Using  the 
average  values  of  four  boundary  points  or  linear /exponential  projection  of 
boundary  points  offers  the  possibility  of  generating  good  starting  iterates. 
Experience  indicates  that  solution  for  both  single  and  multiple  airfoils,  is 
sensitive  to  the  initial  guess  and  distance  to  the  outer  boundary  and  in 
some  cases  a  solution  may  not  be  obtained.  The  solution  often  converges 
slowly. 

5  PARABOLIC  SCHEME  AND  APPLICA¬ 
TION 

5.1  GENERAL  OBSERVATIONS 

In  the  present  study  we  examine  the  feasibility  of  using  parabolic  partial 
differential  equations  for  grid  generation  by  considering  the  pair  of  model 
equations 

HZ*  V)xv  =  HZ,  »?)*«  +  ,  t})Vx{£,  rj)  +  <f(£,  7/)  (57) 

HZ*  v)v*i  =  HZ*  V)t kt  +  c(£>  riVyiZ,  tj)  +  d((,  r?)  (58) 

where  a,b  and  c  can  be  constants  or  some  functions  of  (£,77);  (r,  y)  denote 
the  coordinates  of  the  physical  domain,  (£,77)  the  computational  domain, 
and  (yx,Vy)  source  terms.  The  physical  and  computational  domains  are 
shown  in  fig.  14.  Unlike  the  elliptic  scheme,  in  the  present  approach  the 
inner  boundary  represents  an  initial  condition,  and  the  outer  boundary 
represents  a  constraint  that  effects  the  ./-coordinate  line  distribution.  The 


point  distribution  on  remaining  segments  is  given  and  represents  boundary 
conditions  on  left  and  right  sides  of  the  computational  plane.  Equations 
(57,58)  can  be  discretized  using  backward  differencing  in  the  (timelike)  77- 
coordinate.  Given  the  initial  values  of  x  and  y  at  77  =  0  this  leads  to  a 
tridiagonal  system  to  be  solved  for  each  increment  of  77. 

The  initial  values  are  specified  as 

x(€i  0)  =  ®o(0 

»(C.0)  =  jto(O  (59) 

where  x0(0  and  y0(0  are  the  coordinates  of  the  body  surface.  In  order  to 
examine  the  effect  of  the  source  terms  ,  let  us  set  b  =  d  =  0  so  that 

dx{i,rj) 


drj 


=  -vx  (e,u) 

a 


(60) 


3y((,v)  c,,,, 

— =  (61) 

From  (60,61)  we  see  that  as  77  increases,  the  change  of  x  and  y  is  determined 
by  Vx  and  Vy.  This  implies  that  Vx  and  Vy  should  be  specified  in  such  a 
way  that  x  and  y  change  in  a  desired  direction  and  amount.  The  role  of  x^ 
and  y <#  in  eqs.  (57,58)  may  be  considered  as  smoothing  the  grid  intervals 
in  the  f  direction.  Vx  and  Vy  can  be  approximated  in  (60,61)  by  using 
linear  or  polynomial  interpolation  between  the  inner  and  outer  boundaries. 
The  orthogonality  of  the  grid  lines  may  be  controlled  by  introducing  a 
“fictitious”  outer  boundary,  and  grid  spacing  in  both  the  f  and  77  directions 
can  be  controlled  by  discretization  of  the  governing  equations  (57,58)  on  a 
nonuniform  mesh  in  the  computational  plane. 


5.2  NUMERICAL  DEVELOPMENT  OF 
PARABOLIC  SCHEME 

In  most  coordinate  transformations  used  for  computational  analyses  the 
grid  spacing  in  both  the  £  and  77  directions  on  the  computational  domain 
is  set  to  unity.  However,  this  restriction  is  not  necessary.  When  deriving 
generation  equations,  we  can  assume  that  grid  spacings  are  quasi-  uniform 
on  the  computational  domain.  Once  a  grid  is  generated,  it  may  be  used 
for  flow  calculations  as  if  generated  on  a  uniformly  spaced  grid  on  the 
computational  domain.  Using  Taylor  series  expansion  on  a  nonuniform 
mesh  and  referring  to  fig.  15a,  derivatives  can  be  approximated  as  follows 
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and  the  governing  equations  (57,58)  can  be  discretized  as 
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In  order  to  determine  coefficients  a,b  and  c  and  to  explain  the  proposed 
method  of  controlling  grid  spacing  the  idea  is  to  relate  the  model  equa¬ 
tions  to  some  known  set  of  equations  in  the  computational  plane.  For  this 
purpose,  the  grid  generation  of  elliptic  type  (45),  is  selected.  In  fact,  the 
procedure  that  we  want  to  apply  is  essentially  equivalent  to  deriving  the 
elliptic  difference  equations  (45)  on  the  nonuniform  grids  of  the  computa¬ 
tional  domain  in  which  the  m-th  grid  line  in  the  77-direction  on  a  uniform 
grid  is  moved  toward  the  (m  —  l)-th  grid  line.  Using  the  above  procedure, 
an  arbitrary  number  of  grid  lines  may  be  moved  in  both  the  £  and  77  direc¬ 
tions.  By  denoting  the  distance  ratios  by  F,-  and  gj  in  the  £  and  77  directions, 
respectively,  (fig.  15b),  the  difference  equations  that  yield  the  grid  spacing 
can  be  written  as 


2 A  fS«-i,j  —  Sj,j  ~  si,j  1 

Fi  +  Fi-x1  Fi-i  Fi  J 

org«+l,j+l  -  fr+lj-l  ~  g«-W+l  +  S.-lJ-li  , 
(Fi-i  +  Fi)(gj- 1  +  gj) 

2C  —  Sjj  ^  g»,j-n  ~  si,j  j  _  q 

9i  +  9j- 1  9i- 1  9i 


(64) 


We  then  have  elliptic  PDE’s  (45)  discretized  on  a  quasi-uniform  mesh  in 

(£,*?)• 


Now,  a  marching  grid  generation  equation  scheme  may  be  obtained 
from  (64)  by  replacing  the  coordinates  a\,j+1  and  Vi,j+i  by  known  values 
XOij+i  and  YOij+ 1  obtained  from  linear  interpolation  between  coordinates 
of  the  outer  and  inner  boundaries.  The  ratio  gj  can  be  relaxed  in  practice 
i.e.,  it  may  be  gradually  changed  along  the  grid  line  in  some  functional 
dependence. 

In  the  simplest  case,  when  we  do  not  force  orthogonality  on  the  bound¬ 


aries,  coordinates  XO  and  YO  are  set  to  the  known  values  {XBijmax  , 
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YBi, jmax  )  at  the  outer  boundary.  In  general,  however,  the  marching  grid 
generation  equations  for  both  the  x  and  y  coordinate  may  be  written  in  the 


form: 


2 A  r8i-l,j  Sj,j  .  6j+l,j  ~  8i,j  -I  , 

F  +  F-i  F— i  +  F  J 

^G  r  ^t.J  _  &«,j  i 

9j- 1  +  Gj  g;-i  Gj 
_  r5Q,+i,j-n  —  50,-ij-i  —  g«+i,j-i  +  <s,_  i,y— i  - 
(Fi- 1  4-  F)(0j- 1  +  Gj 

2G  rg»,j-i  SOj,j+i  i 

ffj-i  +  Gj  5j_i  Gj 


where: 


Gj  —  </j  +  Qj+i  +  ••••9  jmax— \  =  tjmax  —  f] j 

9j- 1  =  Vi  “  ty-i  ;  Fi  =  &+1  -  £*  ;  F-i  =  &  -  &-i 

Gj  represents  the  distance  between  grid  lines  (i,j)  and  (*,  JO)  on  the  com¬ 
putational  plane  and  5G,,j  =  XGj,j  or  Y O.  j.  In  order  to  establish  a  relation 
between  discretized  model  parabolic  system  (62,63)  and  system  (65),  let  us 
compare  those  two  systems.  Separating  known  and  unknown  values  we  get 

JL) _ 2si'j  -  (  a  4-  — )s-  |  b  2si-hj  . 

kal  ar’al  +  ar  ka  d  a  o’  AR  (AR  +  AX) 

b  2  si+1J  SB^jmax  _Sx,i- 1 

ahIaFTaZ)  - _c  ao  -aTD~d  (66) 


/  \  _ /  G  G .  2s,- tj  -A  , 

^  F. a  F ;  (F  +  F-i)  v  ^  ^  Gj ;  Gj  +  9j—i  F-!  (F,  +  F-i) 

A  2 Sj+ij  _ _ 2G  j  SOjj+i  8i,j- i  i  _ 

Fj  (F  +  F_i)  Pj-i  +  Gj  Gj  3j-i 
nr^Gi+ij+l  -  50,_ij+j  “  St+W-l  +  Si-1, 

B[  W +  «-,)(»-. +  Gi) - X67) 
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In  both  equations  the  right  sides  are  known  values.  It  can  be  seen  that  they 
have  the  same  form.  So,  we  can  establish  a  relation  between  coefficients  of 
eqs.  (62,63)  and  (65)  as  follows 

b  =  A 


a  =  c  = 


(Gj  +  gj-i) 


»  pr£j2j±W±i  ~  —  -St+ij-i  + 

1  (Fi  +  Fi-MGj  +  gj-x)  J 

Where  the  coefficients  A,B  and  C  are  given  by  (43).  In  fact,  system  (65) 
represents  a  2  x  2  block  tridiagonal  system  that  can  be  evaluated  as  follows: 
defining 

2  A  2  A 

01  ~  +  Fi-i)  5  7  "  Fi(Fi  +  Fi-i) 

0  (1  +  J_) _ — — (_L  +  _L) 

P  Fi  +  Fi-^Fi*  FiJ  Gj  +  9i-iKGj  / 

Dx  =  —B[XOi+ ij+i  —  XOi-ij+i  —  Xi+ i,j_i  4-  x,_iij_i]/[(i?i  -f-  Fi-\)(Gj 

2  c 

+  9j~ i)]  -  r  ,  — [xi,j-i/9j-i  +  XOij+i/Gj] 

Gj  +  9j- 1 


A/  =  Oi+ lj+i  —  y  Oj-lj+l  —  Jfi+ij-l  +  Ifc-lj-lJ/KJ5*-  +  Fi-i)(Gj 

2C 

+  5j-i  )]  -  r  ,  — [y« j-i/ 5j-i  +  ^ A>+i  / Al 
Oj  +  0j-i 

the  governing  equation  in  the  x-direction  becomes 

ax*-ij  +  faij  +  T*i+ij  =  A  (68 

and  similarly  in  the  y- direction: 

+  PVij  +  TVi+ij  =  Dv  (69 
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Equations  (68)  and  ’9)  are  solved  simultaneously  for  all  grid  points  on 
the  j-th  grid  line  by  using  a  tridiagonal  solver  for  each  Xij  and  j h,j-  The 
solution  starts  with  j  =  2,  (next  to  the  inner  boundary)  and  marches  until 
j  =  J  MAX  —  1,  (next  to  the  prescribed  outer  boundary  grid  line).  The 
coefficients  of  A ,  B  and  C  are  calculated  using  the  coordinates  of  adjacent 
grid  lines  that  are  already  generated,  and  the  first  derivatives  that  appear 
in  them  can  be  approximated  by  differencing  as, 


xi-l,j-l 

4  Fi  +  Ft- 1 

(70a) 

-  Vi-u-i 

Fi  +  Fi-r 

(70b) 

XOij+i  —  Zi'j-i 

Si-1  +  Gj 

(71a) 

Y  Oij+i  —  yij-i 

Si-i  +  Gj 

(71b) 

The  non-uniform  grid  spacing  terms  in  the  above  scheme  have  effects  simi¬ 
lar  to  the  spacing  control  terms  in  the  exponential  functions  introduced  by 
Thompson,  Thames  and  Mastin  [3]  for  the  elliptic  scheme.  The  value  of  F 
or  g  and  the  distance  between  the  adjacent  grids  satisfy  approximately  a 
linear  relationship.  Therefore,  if  a  grid  is  generated  by  using  a  known  set 
of  F  and  g  for  all  the  grid  intervals,  and  if  a  different  spacing  distribution 
is  desired  in  the  next  grid  generation  the  values  of  F  and  g  that  satisfy  the 
desired  grid  spacing  can  be  found  by  the  linear  relationship  between  F  and 
g  and  the  grid  spacing  in  the  previous  calculation. 

5.3  GRID  SPACING  CONTROL 

In  computation  of  viscous  flows  or  boundary  layers,  orthogonality  of  grid 
lines  near  the  boundary  surface  is  desirable  to  represent  all  normal  deriva- 


tives  simply  and  accurately.  To  obtain  this,  the  angular  orientation  of  grid 
lines  may  be  controlled  by  using  a  modified  outer  boundary  grid  to  deter¬ 
mine  the  right  hand  side  source  terms.  For  example,  if  the  grid  is  shifted 
in  such  a  way  that  the  new  outer  boundary  points  lie  on  lines  which  are 
perpendicular  to  the  body  surfaces,  the  grid  lines  will  be  away  from  the 
inner  surface  approximately  orthogonally. 

Let  us  first  consider  the  coordinates  ( XOi,3,YOit3 ),  corresponding  to  the 
grid  line  j  =  2.  Referring  to  fig.  16  a  straight  line  A  A  is  extended  outward 
normal  to  the  airfoil  surface  from  the  grid  (i,  1).  A  circular  arc  CC  that 
passes  through  the  point  {XB^jmax  ,YB{tjMAx)  on  the  outer  boundary 
and  has  its  center  at  (i,  1)  is  drawn.  The  coordinates  (XO.,3  and  YOi ,3) 
are  set  to  those  of  the  intersection  S  of  the  two  lines  AA  and  CC. 

In  the  present  work,  the  modified  outer  boundary  is  computed  by  first 
calculating  the  slope  of  the  body  surface  and  the  distance  from  the  body  sur¬ 
face  points  to  corresponding  outer  boundary  points.  The  modified  bound¬ 
ary  is  scaled  so  that  it  lies  at  the  same  distance  from  the  body  surface  as 
the  actual  outer  boundary.  This  is  needed  to  obtain  the  desired  clustering 
or  stretching  of  grid  lines.  Since  the  parabolic  algorithm  generates  the  grid 
lines  progressively  toward  the  outer  boundary  the  modified  grid  is  gradu¬ 
ally  shifted  back  to  the  desired  outer  boundary  point  distribution.  In  this 
manner,  angle  control  is  maintained  at  the  body  surface,  and  a  smooth 
transition  can  be  made  to  the  outer  boundary. 

5.4  EFFECT  OF  THE  COEFFICIENTS  ON  GRID 
GENERATION 
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Here  we  explore  the  possibility  of  generating  grids  with  different  values  of 
coefficients  in  (68,69).  We  consider  the  cases: 
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CASE  1: 


CASE  2: 


CASE  3: 


CASE  4: 


A  =  A(f,ij),  B  =  B((,  rj),  C  =  C((,  rj) 


A  =  A(£,r1),B  =  0,C=C(Z,T)) 


A  =  A(£,ri),B  =  0,C=l 


A  =  0,B  =  B(t,T)),C  =  C(S,r,) 


Case  1  represents  the  method  applied  in  the  present  study.  Coses  2  and 
3  produce  reasonable  grids  (see  fig.  17a  and  fig.  17b,  respectively),  but 
in  these  cases  smoothness  of  the  grid  lines  is  poor.  Case  4  is  interesting 
because  it  provides  a  completely  explicit  grid  generation  method  in  the  rj 
direction.  Grids  from  case  4  are  given  in  fig.  17c  and  we  can  see  that  they 
are  very  similar  to  that  obtained  for  case  1.  Other  combinations  of  coeffi¬ 
cients  may  not  provide  grids  because  boundary  information  is  insufF  ient 
or  due  to  the  changing  type  of  the  governing  equations. 

5.5  PRACTICAL  APPLICATION  AND  RESULTS 

•  SINGLE  AIRFOILS 

In  the  present  work,  all  three  commonly  used  types  of  grids  -  H-type,  O-type 
and  C-type  -  are  generated.  The  complete  procedure  is  given  for  O-type  of 
grids,  and  the  other  two  briefly  explained. 

First  of  all,  let  us  consider  the  transformation  of  a  two-dimensional  dou¬ 
bly  connected  region  D  bounded  by  two  closed  contours  onto  a  rectangular 
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region  D,  as  shown  in  fig.  18  (the  body  contour  and  outer  boundary  are 
transformed,  respectively,  to  the  constant  77  -  lines  forming  the  bottom  and 
top  sides  of  the  transformed  region).  Let  Gi  represent  the  inner  boundary, 
G2  outer  boundary,  D  physical  plane,  and  D\  transformed  plane.  In  order 
to  connect  Gi  and  G2  it  is  necessary  to  make  an  arbitrary  cut  between  G3 
and  G4. 

Conceptually  this  can  be  viewed  as  an  opening  of  the  field  at  the  cut 
and  then  a  deformation  into  a  rectangle.  Now,  let  Gi  map  onto  G\,Gi 
onto  (?2,  G3  onto  G%  and  G4  onto  G\.  G*  and  G£  are  constant  77-lines  in 
the  transformed  plane.  Contours  G 3  and  G4  which  connect  the  contours  G\ 
and  G2  are  coincident  in  the  physical  plane  and  form  left  and  right  sides  of 
the  computational  plane  with  x  and  y  values  given  as  pre-assigned  data. 

To  summarize,  at  this  point  we  are  given  “initial”  conditions  on  the  inner 
boundary,  a  constraint  condition  for  the  outer  boundary  and  boundary 
condition  on  cut  as  follows 
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(72) 

(73) 

(74) 


The  functions  /i(£>  *7),  /2(£,  v)i  /3(C)  77),  /4(f)  tj)  are  determined  using  the 
known  shape  of  the  inner  and  outer  boundaries  Gi  and  G2,  and  the  spec¬ 
ified  distribution  of  £  thereon,  while  the  function  /5  is  determined  by  the 
shape  of  the  cut,  and  distribution  of  £  lines  thereon. 

The  boundary  fitted  coordinate  system  generated  by  solving  (68,  69) 
has  a  constant  77  line  coincident  with  each  boundary  in  the  physical  plane. 
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The  £  =  constant  lines  may  be  specified  as  desired  around  boundaries, 
since  the  assignment  of  the  f  -  values  to  the  ( x,y )  boundary  points  via 
the  functions  /i,/2,/3  and  f4  is  arbitrary.  Control  of  the  radial  spacing  of 
the  £  =  constant  lines  is  very  important  because  of  its  influence  on  non- 
uniform  grid  spacing  terms  in  the  finite  difference  equations  (65).  Knowing 
boundary  conditions  (73,  74)  and  initial  conditions  (72)  we  can  calculate 
the  grid  by  marching  from  inner  to  outer  boundary. 

Figures  19  and  20  show  O-type  grids  about  NACA  0012  and  Karman- 
Trefftz  airfoils  generated  in  this  way.  The  grid  points  on  the  outer  boundary 
are  equally  spaced  along  a  circle  with  radius  of  1.5  x  chord  length.  Orthog¬ 
onality  control  is  applied  in  the  vicinity  of  the  airfoil  surface.  Grid  size  for 
NACA  0012  is  61  x  28  (CPU  time  0.258  sec.)  and  for  Karman  -  Trefftz 
airfoil  41  x  16  (CPU  time  0.189  sec.).  Fig.  21  shows  an  O-grid  generated 
about  a  NACA  0012  airfoil  with  strong  clustering  in  vicinity  of  the  air¬ 
foil  surface.  This  type  of  coordinate  line  clustering  can  be  considered  for 
viscous  and  layer  problems. 

Figure  22  shows  H-type  grids  without  orthogonality  at  the  airfoil  surface 
and  generated  in  the  same  way.  The  horizontal  line  of  symmetry  is  used  as 
initial  data.  Distribution  of  points  on  the  outer  boundary  is  the  same  as  on 
the  line  of  symmetry.  Essentially,  the  grid  is  equivalent  to  that  obtained  by 
an  algebraic  grid  generation  scheme  applied  by  Johnson  [34]  for  transonic 
flow  calculation  about  a  single  airfoil  in  channel  conditions.  The  grid  size 
is  41  x  20  and  solution  is  obtained  after  0.126  sec  of  CPU  time.  The  grids 
above  the  airfoil  were  first  generated  starting  from  the  airfoil  surface  to 
the  top  boundary,  and  then  grids  below  the  airfoil  were  generated  by  the 
same  procedure.  The  best  grid  resolution  was  obtained  by  specifying  a 
nonuniform  spacing  function  F  in  the  computational  plane  corresponding 
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to  Aar  in  the  physical  plane.  Fig.  23  shows  H-grids  about  a  NACA  0012 
airfoil  with  strong  clustering  and  orthogonality  of  grid  lines  in  the  vicinity 
of  the  airfoil  surface.  Obviously,  both  grids  can  be  used  for  flow  analysis 
about  an  airfoil  in  wind-tunnel  conditions. 

Fig.  24  indicates  C-type  grids  similarly  generated.  In  this  case  the  wake 
position  was  taken  to  define  the  initial  condition  for  resolution  of  grid  lines 
from  the  trailing  edge  of  the  airfoil  in  the  downstream  direction.  The  outer 
boundary  is  uniformly  discretized  and  the  point  distribution  downstream 
from  the  trailing  edge  is  identical  to  that  on  a  wake.  The  grid  size  is  73  x 
16  and  CPU  time  0.231  sec. 

Opening  the  field  at  the  cut  (fig.  25)  the  two  members  of  the  pair  of 
segments  forming  the  branch  cut  are  similarly  directed  in  the  transformed 
region,  and  consequently  points  located  at  a  vertical  distance  below  the 
segment  1-2,  at  a  horizontal  distance  to  the  left  of  point  2,  coincide  with 
points  at  the  same  vertical  distance  above  the  segment  4-3,  at  the  same 
horizontal  distance  to  the  left  of  point  3.  In  this  case,  £  varies  to  the 
right  on  the  upper  side  of  the  cut,  but  to  the  left  on  the  lower  side.  The 
direction  of  variation  of  77  also  reverses  at  the  cut,  so  that  although  the 
type  and  shape  of  both  lines  are  continuous  across  the  cut,  the  direction  of 
variation  reverses  there. 

In  all  cases  shown,  for  O-type  and  C-type  grids  the  nonuniform  spacing 
functions  F  were  set  to  ds  =  \J Ax'-  +  A yf  on  the  j  —  1  grid  line,  while 
nonuniform  spacing  of  the  77  coordinate  was  enforced  by  use  of  trigonometric 
sine  functions  to  produce  clustering  in  desired  regions. 
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The  basic  ideas  and  procedure  introduced  for  single  airfoils  can  be  extended 
to  regions  containing  more  than  one  body,  i.e.,  multiconnected  or  multi  - 
body  configurations.  An  example  of  the  transformation  for  two  airfoils  is 
given  in  figure  26.  The  bodies  are  connected  with  one  arbitrary  cut  and 
an  additional  arbitrary  cut  joining  one  of  the  body  contours  to  the  outer 
boundary.  The  physical  plane  contours  G\  —  G&  map  respectively  onto  the 
contours  G[  —  in  the  transformed  plane.  The  conceptual  opening  here  is 
as  follows:  the  pairs  of  segments  (1-2, 7-8)  and  (3-4, 5-6)  are  the  branch  cuts, 
which  form  re-entrant  boundaries  in  the  transformed  plane.  In  this  case, 
points  outside  the  right  side  of  the  transformed  region  coincide  with  points 
inside  the  left  side,  and  vice  versa.  The  coordinate  type  and  direction  are 
continuous  across  the  cut.  Points  below  the  bottom  segment  3-4,  to  the 
left  of  point  4,  coincide  with  points  above  the  segment  5-6  to  the  right  of 
point  5  in  the  transformed  plane.  As  those  arbitrary  cuts  represent  pre¬ 
assigned  data  for  grid  generation  in  the  physical  plane,  there  are  a  number 
of  other  possibilities  for  placement  of  the  two  cuts  on  the  boundary  of  the 
transformed  region. 

In  the  present  work,  the  procedure  was  applied  to  two  NACA  0012 
airfoils  with  flap  at  different  angles.  Examples  with  flap  angle  at  0  degrees 
and  25  degrees  are  shown  in  figures  27  and  28,  respectively.  Figure  29  gives 
more  detailed  grid  line  distribution  for  two  NACA  0012  airfoils  with  flap 
angle  of  25  degrees.  The  grid  size  is  87  x  30  and  CPU  time  0.623  sec.  The 
results  show  very  good  resolution  of  grid  lines  and  clustering  effects  of  the 
nonuniform  spacing  terms. 


•  APPLICATION  TO  ARBITRARY  CONFIGURATIONS 


The  method  is,  with  the  same  basic  ideas,  applied  further  to  cascade  grids 
and  configurations  with  two  circles  and  a  circle  with  an  airfoil. 

The  first  case  demonstrates  the  ability  of  the  parabolic  method  to  satisfy 
outer  boundary  conditions  a  very  small  distance  from  the  inner  boundary. 
In  this  case  the  outer  boundary  was  formed  midway  between  neighboring 
airfoils.  In  fact,  this  cascade  grid  represents  a  composite  system  that  can  be 
extended  upstream  and  downstream  to  infinity  by  two  independent  Carte¬ 
sian  systems  similar  to  that  given  by  Eiseman  [7].  The  main  advantages  of 
the  composite  system  are  that  it  can  be  used  for  highly  -  staggered,  closely 
-  spaced  airfoils  and  avoids  severe  mesh  distortion  or  growth  that  would 
occur  at  higher  upstream  or  downstream  extensions  of  the  basic  O-grids. 
The  mesh  in  fig.  30a  covers  everything  except  triangular  regions  at  the 
corner  points  that  result  from  the  demand  for  nonsingularity.  Obviously, 
this  system  is  to  be  applied  only  when  these  regions  are  sufficiently  small 
and/or  are  located  in  places  where  the  solution  is  slowly  varying.  An  al¬ 
ternative  approach  is  to  relax  the  demand  for  nonsingularity  and  fill  in  the 
uncovered  regions  as  in  fig.  30b. 

The  configuration  with  two  circles  (fig.  31)  shows  the  clustering  ability 
of  the  method.  Essentially  the  grids  of  this  problem  have  very  similar 
behavior  to  the  example  using  an  elliptic  generator  given  by  Thompson, 
Thames  and  Mastin  [16]. 

The  configuration  containing  circle  and  airfoil  shapes  (fig.  32)  indicate 
the  ability  of  the  parabolic  scheme  to  handle  grid  generation  in  a  field  with 
multiple  bodies  having  very  different  geometry.  The  example  can  be  applied 
for  grid  generation  about  an  arrow  wing  -  body  configuration. 
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5.6  COMPUTATIONAL  EFFICIENCY 

The  parabolic  scheme  produces  viable  grids  in  one  outward  “iteration” 
sweep.  It  uses  a  tridiagonal  solver  that  is  very  time  efficient.  In  order 
to  investigate  the  computational  efficiency  of  the  parabolic  scheme  as  a 
function  of  grid  size  we  compare  results  with  those  obtained  by  the  elliptic 
scheme  with  clustering  functions.  Results  are  given  in  Table  1  and  show 
that  the  parabolic  scheme  is  highly  efficient.  Morever  it  exhibits  very  little 
dependence  on  the  distance  of  the  outer  boundary.  In  all  cases  for  a  single 
NACA  0012  airfoil  CPU  time  includes  forming  the  inner  boundary,  outer 
boundary  and  grid  generation.  In  all  other  cases  CPU  time  is  for  forming 
of  the  outer  boundary  and  grid  generation.  All  presented  examples  were 
run  on  the  CDC  170/750  Dual  Cyber  Computer  of  the  University  Texas  at 
Austin. 

5.7  LIMITATIONS  OF  PARABOLIC  SCHEME 

We  find  that  the  most  important  and  the  most  sensitive  factors  influencing 
the  method  are  the  distribution  of  points  on  cuts  and  specification  of  the 
nonuniform  spacing  terms  Fi  and  gr  In  order  to  obtain  good  resolution  of 
grid  lines  in  the  cut  region  for  varying  nonuniform  spacing  terms  Fi  and  gj 
in  (68,  69)  with  j,  the  Fi  and  gj  should  be  given  as  functions  of  the  point 
distribution.  In  this  case  clustering  of  the  grid  lines  in  the  entire  physical 
field  is  controlled  by  grid  point  distribution  on  the  cuts.  Otherwise  we  can 
have  excessive  skewness  of  the  grids  near  the  cuts.  If  we  use  constant  values 
of  Fi  and  gj ,  to  obtain  different  clustering  of  the  grid  lines  we  need  several 
sweeps  through  the  computational  plane  changing  Fi  and  gj  in  each  sweep. 

In  application  to  multiple  bodies  we  have  additional  cuts  that  connect 


|g 

m 


m 

M 

mm 

M 

mm 

vvv' 

m 


/Am 


f:#! 


48 


the  bodies.  In  the  case  of  asymmetric  bodies  we  found  that  uniform  spacing 
of  points  on  cuts  is  most  appropriate,  since  otherwise  overlapping  may 
occur.  For  bodies  symmetrically  situated  in  the  physical  plane,  these  cuts 
are  not  critical  and  the  point  distribution  on  them  may  be  arbitrary. 

5.8  CLUSTERING  OF  POINTS  IN  GIVEN 
REGIONS 

The  distribution  and  clustering  of  points  along  the  airfoil  contour,  cuts, wake 
and  “inlet”  and  “outlet”  of  H-type  grids  are  controlled  using  a  trigonometric 
sine  function.  We  have  two  known  parameters;  the  coordinate  location  of 
a  point  and  the  number  of  grid  points  that  will  be  located  up  to  the  known 
point.  The  selected  function  is 

y  =  L[x  +  FAC(Nirx)]  (75) 

where  x  linearly  varies  from  zero  to  the  length  of  the  interval,  FAC  is  a 
factor  that  determines  the  amount  of  clustering,  N  is  1  or  2  (for  N  =  1 
we  have  a  single  and  for  N  =  2  a  double  sine  wave),  and  L  is  the  physical 
length  of  the  interval. 

6  CONCLUSIONS 

The  present  work  shows  that  high  quality  grids  can  be  generated  by  a  mar¬ 
ching  solution  of  parabolic  partial  differential  equations.  The  parabolic  dif¬ 
ferencing  scheme  has  been  implemented  for  the  generation  of  two-dimensional 
grids  about  single  and  multiple  bodies.  The  method  has  demonstrated  good 
spacing  and  angle  control,  and  ability  to  generate  smooth  grids  for  finite 
difference  computations  in  CFD. 


The  approximate  orthogonality  of  the  grids  is  enforced  by  introducing 
a  fictitious  outer  grid  contour.  This  shifted  outer  grid  contour  is  computed 
by  first  computing  the  the  slope  of  the  airfoil  surface  and  the  distance  from 
inner  boundary  points  to  corresponding  outer  boundary  points. 

On  the  basis  of  comparison  of  results  with  those  obtained  by  an  elliptic 
scheme,  the  method  is  demonstrated  to  have  high  computational  efficiency 
and  gives  good  resolution  of  grid  lines.  The  present  method  appears  stable. 
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Part  II 

POTENTIAL  AIRFOIL 
ANALYSIS  AND  DESIGN 


7  INTRODUCTION 

The  most  popular  and  practical  applicable  methods  for  potential  airfoil 
analysis  and  design  are  so  called  surface  singularity  or  panel  methods.  The 
basic  ideas  of  the  panel  methods  introduced  in  the  solution  of  arbitrary 
potential  flow  problems  involve  combining  classical  potential  theory  with 
contemporary  numerical  techniques.  The  classical  theory  provides  a  means 
to  reduce  flow  problems  to  a  surface  integral  equation  relating  boundary 
conditions  to  an  unknown  singularity  distribution.  Numerical  techniques 
are  then  used  to  calculate  an  approximate  solution  of  the  integral  equa¬ 
tions.  The  procedure  involves  representing  flow  boundaries  by  surface  el¬ 
ements  (panels)  on  which  potential  flow  singularities  are  distributed.  The 
physical  boundary  condition  is  that  the  normal  velocity  is  specified  on  all 
surfaces  so  that  the  mathematical  formulation  becomes  a  Neuman  problem 
for  Laplace’s  equation. 

For  a  given  geometry  of  an  airfoil  section  the  potential  flow  analysis 
methods  provide  surface  velocity  and  pressure  distribution  and  thus  the  lift, 
the  pitching  moment  and  other  aerodynamic  characteristics.  The  methods 
should  be  reliable  and  accurate  because  they  are  used  as  the  first  step  in  the 
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design  of  new  wing  sections.  Generally,  exact  (i.e.,  conformal  mapping)  and 
surface  singularity  methods  provide  a  good  means  for  single  airfoil  poten¬ 
tial  flow  analysis  and  design.  However,  for  multicomponent  airfoils  exact 
methods  produce  mainly  test  cases,  and  the  only  practical  design  methods 
are  based  on  surface  singularity  techniques.  In  a  design  mode  we  usually 
specify  velocity  or  pressure  distribution.  Therefore,  the  performance  of  the 
system  specified  by  design  or  inverse  solution  gives  the  geometry  that  will 
produce  that  performance.  The  inverse  design  is  basically  an  iterative  pro¬ 
cess  that  provides  adjustments  to  both  designed  geometry  and  the  required 
surface  velocities  at  each  iteration.  So,  we  can  include  different  constraints 
that  frequently  occur  in  the  design  of  airfoil  sections. 

7.1  REVIEW  OF  POTENTIAL  METHODS 

The  most  widely  used  panel  method  is  that  of  Hess  and  Smith  [35]  based 
on  distribution  of  sources  and  sinks  on  the  airfoil  surface  combined  with 
a  vorticity  distribution  to  generate  circulation.  Improved  solution  based 
on  a  higher-order  panel  formulation  was  developed  by  Hess  [36,  37,  38]. 
The  main  idea  in  development  of  higher  order  panel  methods  was  to  reduce 
panel  density  for  a  given  solution  accuracy  and  to  reduce  computing  cost. 
However  these  advantages  have  not  been  achieved.  For  example,  Maskew 
[39]  demonstrates  that  a  low-order  panel  method  based  on  piecewise  con¬ 
stant  doublet  and  source  singularity  applied  to  general  configurations  for 
comparable  density  of  control  points  give  comparable  accuracy  to  higher- 
order  solutions.  Bristow  and  Grose[40]  introduce  panel  methods  based  on 
Green’s  third  identity  with  constant  or  linear  source  distribution  and  linear 
vorticity  distribution  on  flat  panels  and  demonstrate  very  efficient  inverse 
and  mixed  analysis/design  for  single  and  multicomponent  airfoils.  For  two- 
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dimensional  configurations  a  very  popular  method  is  the  stream  function 
approach  based  on  constant  or  linear  distribution  of  vorticity.  For  exam¬ 
ple,  Ormsbee  and  Chen  [41]  use  this  method  for  multi-element  airfoil  design 
and  use  third-order  Langrangian  interpolation  to  determine  surface  points. 
Further,  Kennedy  and  Marsden  [42]  introduce  an  additional  control  point 
a  finite  distance  from  the  trailing  edge  to  satisfy  the  Kutta  condition  and 
determine  surface  points  using  cubic  spline  or  linear  projection.  This  type 
of  Kutta  condition  provides  reduction  in  panel  density  and  gives  very  good 
results.  The  method  is  also  used  as  a  starting  solution  in  design  of  single 
airfoils  in  viscous  incompressible  flow  by  Dutt  and  Sreekanth  [43]  and  for 
analysis  and  design  of  single  airfoils  in  transonic  flow  by  Greff  and  Man¬ 
tel  [44] 

Instead  of  superposition  of  potentials  due  to  surface  singularities,  a  dif¬ 
ferent  type  of  boundary  integral  equation  can  be  obtained  by  applying  an 
appropriate  Green’s  formula  and  fundamental  solution  to  relate  integrals 
over  the  interior  of  the  domain  to  integrals  on  the  boundary.  Then,  in¬ 
troducing  a  finite  element  expansion  on  a  discretization  of  the  boundary 
domain,  an  approximate  solution  of  the  boundary  integral  equations  can  be 
obtained.  This  procedure  has  been  termed  the  boundary  element  method 
(Brebbia  [45]).  The  method  has  been,  with  theoretical  development,  suc¬ 
cessfully  applied  to  lifting  airfoil  calculations  by  Carey  and  Kim  [46].  The 
major  conceptual  difference  between  panel  and  boundary  element  meth¬ 
ods  is;  the  panel  methods  are  based  on  superposition  of  surface  singularities 
with  discrete  satisfaction  of  boundary  flux  conditions  while  in  the  bound¬ 
ary  element  methods  one  uses  a  finite  element  expansion  and  a  discrete 
approximation  of  the  boundary  integral  equations. 

Conformal  mapping  techniques  have  been  successfully  applied  by  Halsey 
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[9]  for  analysis  of  multi-element  airfoils.  Any  number  of  airfoils  are  trans¬ 
formed  to  the  same  number  of  circles  by  successive  application  of  a  method 
for  mapping  a  single  body  to  a  unit  circle.  Then,  the  flow  about  multiple 
circles  is  analysed  using  multiply-reflected  doublets  and  vortices.  Saddhoo 
and  Hall  [47]  use  the  method  of  images  and  calculate  the  analytical  solution 
for  an  inviscid  incompressible  flow  past  four  arbitrary  circles  which  are  then 
conformally  mapped  onto  airfoil  sections. 

In  the  present  work  we  use  the  stream  function  approach  of  Kennedy 
and  Marsden  [42].  However,  in  the  design  mode  we  calculate  the  actual 
position  of  the  trailing  edge  and  then  use  linear  interpolation  to  determine 
surface  points.  In  this  way  it  is  possible  to  avoid  saw-shaped  airfoils  that 
can  result  from  propagation  of  small  errors  due  to  selection  of  additional 
trailing  point  as  the  actual  trailing  edge  of  the  airfoil.  The  method  is 
compared  with  several  panel  methods  in  analysis  mode  and  demonstrated 
approximately  to  be  of  second  order  accuracy. 

7.2  GENERAL  CONSIDERATIONS 

The  first  step  in  the  solution  of  potential  flow  over  single  or  multi-element 
airfoils  is  to  define  elements  that  describe  the  airfoil  surface.  We  select  a 
certain  number  of  points  on  the  airfoil  surface  and  connect  those  points 
with  straight  or  curved  lines  which  define  the  panels  (fig.  33). 

Next,  w e  (1)  represent  the  body  surface  and  its  wake  by  a  distribution  of 
sources,  doublets,  and/or  vortices  of  unknown  strength,  (2)  parameterize 
the  singularity  strength  (i.e.,  represent  it  by  a  polynomial),  (3)  enforce 
an  appropriate  boundary  condition  for  the  velocity  field  at  control  points, 
and  (4)  solve  the  resulting  system  of  linear  algebraic  equations  in  unknown 
singularity  strength. 
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In  order  to  be  able  to  treat  lifting  bodies  it  is  necessary  to  introduce 
circulation,  the  strength  of  which  is  determined  by  the  Kutta  condition. 
Using  an  analysis  of  the  potential  flow  near  a  sharp  convex  edge  (Moran 
[48]),  there  are  two  possibilities:  (a)  the  velocity  at  the  trailing  edge  is 
infinite;  and  ( b )  the  flow  leaves  the  sharp  trailing  edge  of  an  airfoil  section 
in  a  smooth  fashion.  From  the  requirement  of  finite  velocity  at  the  trailing 
edge,  the  second  alternative  gives  the  following  practical  applications  of 
the  Kutta  conditions:  (1)  the  streamline  that  leaves  a  sharp  trailing  edge 
is  an  extension  of  the  bisector  of  the  trailing  edge,  (2)  the  flow  speeds  on 
the  upper  and  lower  surfaces  near  the  trailing  edge  are  equal  at  equal  dis¬ 
tance  from  the  trailing  edge,  and  (3)  if  the  trailing  edge  is  not  cusped,  the 
flow  stagnates  there.  The  most  commonly  used  Kutta  condition  in  surface 
singularity  methods  is  that  airflow  leaves  the  trailing  edge  smoothly.  This 
Kutta  condition  can  be  obtained  by  equating  the  trailing  velocities  over  the 
upper  and  lower  surface  elements  adjacent  to  the  trailing  edge, 

vh  +  vtN  =  0  (76) 

This  form  of  Kutta  condition  results  in  zero  loading  on  elements  nearest 
the  trailing  edge.  In  order  to  minimize  the  error  due  to  this,  we  have  to 
use  very  small  elements  in  the  trailing  edge  region.  Other  possibilities  are 
to  use  an  additional  control  point  a  finite  distance  from  the  trailing  edge  or 
to  approximate  the  tangential  velocity  component  at  the  trailing  edge  by 
two-  or  three-point  quadratic  extrapolation. 
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8  MATHEMATICAL  DEVELOPMENT 
OF  PANEL  METHODS 


Consider  compressible,  steady  and  irrotational  flow  of  inviscid  fluid  in  three 
dimensions.  The  fact  that  the  flow  is  irrotational  implies  the  existence  of 
a  velocity  potential  V$  =  V.  The  continuity  equation  then  reduces  to  the 
full  potential  equation 
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For  slightly  compressible  flow  we  separate  terms  to 


(77) 


A$  =  /($) 


(78) 


Where  /  corresponds  to  the  nonlinear  terms.  A  simple  Taylor  iterative 
method  then  involves  repetitive  solution  of  the  Poisson  equations 

=  /(«(*-*))  (79) 


Further,  for  incompressible  flow  we  have  the  familiar  linear  problem 


A$  =  0 


(80) 


or  in  terms  of  stream  function  $ 


A«  = 


a2$  a2# 
ax2  +  ay2 


Once  $  or  $  is  determined,  the  velocity  field  can  be  found  as 
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d$  d<u 


(82b) 
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In  order  to  solve  a  potential  flow  problem  for  incompressible  flow,  eq.  (80) 

or  (81)  should  be  discretized  with  specified  boundary  conditions.  We  first 
convert  the  governing  equation  into  an  integral  form.  It  is  known  that  any 
incompressible  flow  can  be  represented  by  a  distribution  of  sources  and 
vortices  over  its  boundary  surface.  To  make  the  potential  $  single  valued, 
we  employ  a  branch  cut  to  the  far  field  (infinity).  The  body  surface,  outer 
boundary  and  cut  surface  are  shown  in  fig.  34. 

The  associated  boundary  conditions  are 


£-0 

on 


on  the  airfoil  surface  Sb  and 


V$  —*  Voo  as  r2  =  x1  +  y2  — > ►  oo  (84) 

in  the  far  field.  Voo  is  the  specified  uniform  flow  at  infinity.  On  the  outer 
boundary  So,  the  boundary  condition  (84)  can  be  approximated  as 

k'v-*  (85) 

Where  n  is  the  unit  vector  outward  and  normal  to  Sb- 

For  lifting  airfoil  problems,  the  unknown  circulation  T  is  determined 
from  the  Kutta  condition.  For  any  simple  closed  contour  C  enclosing  the 
airfoil  the  circulation  T  is  defined  as 

r  =  £  u  ds  =  jf  V$ds  =  $(U)  -  $(W)  (86) 

Where  u  is  the  velocity  and  U ,  W  are  a  pair  of  adjacent  points  on  either 
side  of  branch  cut  Sc-  Across  the  cut  we  have 

[3>]  =  T  on  Sc  (87) 
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and  the  velocity  is  continuous,  so 


d$  n  c 

-t—  =0  on  be 

on 


Where  [  ]  denotes  the  jump  across  the  branch  cut.  Now,  let  us  introduce  a 
source  of  unit  strength  at  some  point  P  in  the  integration  domain  A  (i.e., 
flow  field).  In  this  case  the  fundamental  solution  of  equation  (80)  satisfies 
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A  x  =  S{X-t,Y-ri)  (89) 

and  is  x  —  jvr/nr,  i.e.,  the  potential  of  a  source  of  unit  strength.  Where 
r  =  J(X  -  f)2  +  (Y-  t,)2 
6  is  the  Dirac  delta  function,  and 

£,t]  are  local  coordinates  of  an  arbitrary  point  in  the  integration 
domain  A. 

Multiplying  eqs.  (78)  and  (89)  by  x  and  $,  respectively,  and  subtracting, 
we  get 

$V2x  —  xV2$  =  6$  —  xf  (90) 

The  divergence  of  the  left  hand  side  of  eq.  (90)  is 
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V4>  Vx  +  $  V  Vx  —  xV  V$  —  Vx  V$  =  (5$  —  xf 


V  ($Vx  -xV$)  =  6$  -  xf  (92 

Integrating  over  the  area  A  and  applying  the  divergence  theorem,  we  get 
Ja{6$  -Xf)dA  =  -  Jsn{$VX-xV$)dS  (93 
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Integrating  the  right-hand  side  of  eq.  (93)  by  parts 

L6*dA=Js{xj;-*t)ds+hxfdA  (94> 

and  evaluating  the  right  side  of  eq.  (94),  we  get  the  basic  equation  for  panel 
formulation  as  follows 

w*-') -/,<*£ (95) 

where  the  surface  S  consists  of  three  components  S  =  Sb  +  Sc  +  So 

Sb  is  the  body  surface  immersed  in  the  flow 

Sc  is  a  two-sided  surface  on  the  cut 

So  is  the  surface  of  the  outer  boundary 

and  0  =  ■—  (9  is  the  interior  angle). 

Hence,  the  governing  equation  becomes 

mi,ri)=L  „  Ax^-^)dS  (96) 

JS&  4*  Sc  4*  So  OTl  C/Tl 

Where,  0  =  0  if  (£,  rf)  is  outside  of  the  boundary,  0  =  |  if  (f ,  rj)  is  on  a 
smooth  boundary  and  0  =  1  if  (C,t?)  is  inside  the  boundary. 

8.1  POTENTIAL  BASED  PANEL  METHOD 

•  TOTAL  POTENTIAL 

If  $  in  eq.  (80)  represents  a  total  potential,  then  on  the  outer  boundary 
So  we  have 

$  =  VooX  cos  a  +  VooY  sina  (97) 


where  a  is  the  angle  of  incidence  of  the  flow. 


Hence,  the  line  integral  becomes 


(3$  =  Voo(X  cosa +  Y  sina)  +  [  (98) 

JsB+sc  on  on 

Since  d$/dn  =  0  on  a  solid  surface,  integrating  by  parts,  we  get 

0$  =  Voo{X  cos  a  +  Y  sin  a)  +  [  Q^rdS  -  [  (99) 

JsB  On  Jsc  on  v  ' 

where  [  $  ]  =  $+  —  represents  the  potential  jump  across  the  wake  that  is 
equal  to  the  circulation  about  the  airfoil.  In  fact,  may  be  considered  as 
a  doublet  of  unknown  strength  so  that  equation  (99)  represents  a  doublet 
panel  integral  equation  that  is  important  for  the  development  of  vortex 
panel  methods  used  here  in  numerical  experiments. 

•  PERTURBATION  POTENTIAL  METHOD 

If  $  in  eq.  (80)  represents  the  perturbation  potential,  then  velocity  can  be 


written  as 


V  =  Uoo  + 


(100) 


In  this  case  the  potential  vanishes  on  So,  so  that  the  line  integral  becomes 


0*  =  J  (*jL-^x)dS=  [  rx) 

JSb+Sc  On  On  JSg+Sr  On 


f  (*£ 

sB+Sc  on 


(101) 


Where  cr  =  represents  source  strength  and  is  given  as  a  boundary  con¬ 


dition  on  the  airfoil  surface  as 


d*  -t7 
— —  =  —  n  Voc 
on 

Equation  (101)  represents  a  source-doublet  panel  equation. 


VORTEX  METHODS 


(102) 


w '  *< 
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A  surface  doublet  distribution  of  density  p,  can  be  replaced  by  an  equivalent 
surface  vortex  distribution.  So,  vortex  based  panel  methods  can  be  devel¬ 
oped  using  eq.(94).  Integrating  by  parts  and  setting  7  =  | j,f  | *dS  =  0 
and  U  =  0  (since  [  $  ]  is  constant  along  the  wake),  we  get 

0${£,'n)  =  V^X  cosa +  Y  sina)  -  f  y0  dS  (103) 

J  Sb 

or  in  terms  of  a  normal  velocity  at  each  control  point  and  with  /3  =  the 
vortex  integral  equation  becomes 

^  t  _ 


=  -2L7Sa^c 


(104) 


Using  the  relation  between  doublet  and  vortex  distributions  we  can  develop 
the  basic  equation  for  source- vortex  panel  methods.  Recalling  the  doublet 
and  source  equations  (99,  101)  and  knowing  the  potential  jump  [  $  ]  = 
$+  —  and  =  0  on  the  wake,  we  get 


(106) 


*  “  L  *lTndS  +  IsJ  *  'TndS  ~  L  I**™  (105) 

Integrating  by  parts, 

/?$  =  —/  jffds  —  f  <?xdS  (106) 

JSb  JSg 

where  a  —  is  the  unknown  source  strength  and  7  =  |^  the  unknown 
vortex  strength. 

As  a  function  of  the  normal  velocity  at  each  control  point, 

-  r  00  r  dx 

T  r  a  I  v/*'  jp  a  /  v  A  » /-»  /i 


nVoo  =  —2  f  7j?-dS  -  2  /  A 
JsB  on  J sc  on 


(107) 


To  determine  the  appropriate  boundary  integral  solution  to  (80)  or  (81)  it 
is  necessary  to  discretize  the  basic  panel  equations  and  introduce  the  Kutta 
condition  as  a  constraint.  A  detailed  discretization  for  the  panel  method 
with  linear  distribution  of  vorticity  and  a  curved  panel  method  is  given  in 
appendices  A  and  B ,  respectively. 
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9  STREAM  FUNCTION  APPROACH 


As  there  is  no  normal  velocity  at  a  solid  surface,  each  solid  surface  is  a 
streamline  of  the  flow  and  on  the  corresponding  streamline  $  is  a  constant. 
Thus,  for  a  multi-airfoil  problem  the  boundary  conditions  for  eq.  (81)  on 
an  airfoil  component  k  can  be  written  as 

®  =  A*  (108) 

For  computation  it  is  convenient  to  have  all  equations  in  a  nondimensional 
form,  where  distances  are  dimensionless  with  respect  to  the  chord  length 
c,  velocities  with  respect  to  the  free  stream  velocity  14o,  and  the  stream 
functions  with  respect  to  the  product  c.  For  a  uniform  flow  incident  to 
the  positive  x  -  axis  at  an  angle  of  attack  a  the  dimensionless  form  of  the 
stream  function  becomes 

V’/fc  =  ysin  a  —  xcos  a  (109) 

When  the  airfoil  surface  is  replaced  by  a  vortex  sheet,  the  sum  of  the  stream 
function  for  a  uniform  stream  and  the  stream  function  for  the  vortex  sheet 
should  be  constant  on  the  airfoil  surface.  This  can  be  represented  by  the 
integral  equation  (fig.  35)  as 

r/'jt  =  y{S)coso:  —  x(S)sina  —  J  j(S  )ln(S,S  )ds  (110) 

where: 

is  the  unknown  stream  function  value  on  k  -  th  airfoil  section 
7(5')  is  the  vorticity  strength  at  arbitrary  point  5 , 
r(5,  S')  is  the  distance  between  points  5  and  5 
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x(S),y(S)  are  coordinates  of  the  point  of  interest  S 

s,s'  curvilinear  coordinates  measured  along  the  airfoil  surface  starting  at 
the  trailing  edge 

To  solve  equation  (110),  the  airfoil  surface  is  divided  into  N  small  sur¬ 
face  elements  -  panels,  and  the  integral  is  approximated  by  a  summation. 
Applying  eq.  (110)  at  a  control  point  C,  we  obtain 
N  x  , 

A  +  Y'  —  /  y(Sj)lnr(Ci,  Sj)dsj  =  yacos  a  -  xCisina  (111) 

jm  1  **  J‘ 

Assuming  that  we  have  N  control  points  Ci,  the  problem  of  potential  flow 
over  an  airfoil  section  is  reduced  to  that  of  solving  these  N  simultaneous 
equations. 

The  most  immediately  required  result  is  the  velocity  distribution  on 
the  airfoil  surface.  Since  the  velocity  inside  the  airfoil  is  equal  zero,  the 
discontinuity  in  tangential  velocity  across  a  vortex  sheet  is  equal  to  the 
density  of  the  vortex  sheet.  This  implies  that,  in  solving  eq.  (Ill)  for 
7 (Sj),  we  directly  obtain  the  velocity  distribution  on  the  airfoil  surface. 

In  application  of  this  method,  we  assume  that  we  have  straight  panels 
and  constant  (or  linear)  distribution  of  7 j  over  each  element.  Then  we  get 


A  +  =  RHSi  ( i  =  l,2....iV) 

i- 1 


(112) 


where  is  the  influence  coefficient  of  the  element  j  on  control  point  i  and 
RHSi  the  right  hand  side  of  equation  (111)  evaluated  at  control  point  i. 
Using  the  notation  of  fig.  36,  the  influence  coefficients  can  be  written 

r* 7+1 

r  ’  '  'SuS^ds 
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/C,j  =  ^[-fci/n(rx)  +  b2ln{r2)]  -  ^ 

+  ^.[tan-1^- -  *an-1^-]  for  i^j 

-2?r  1)3  03 


Kij  =  ^[(n(^  -  !]  f°r  *  =  3 


(114) 


(115) 


where 


Sj 

r\  =  (xj  -  xCi )2  +  (yj  -  VCi)2 

T-l  -  (Xj+l  -  XCif  +  {Vj+l  -  VCi? 

61  =  “  xc^xi+l  -  *i)  +  (w  “  VCi)(Vj+ 1  -  »)] 

62  =  X^Xj+1  ”  *<*)(*i+i  ~  xi)  +  (»+!  “  Jto)(%+1  “  Vj)] 

63  =  “  Xc-^^+1  ~  w)  “  (»  ”  yci)(xj+i  -  Xj)] 

The  vortex  strength  7;-  at  the  intersection  of  two  panels  is  determined  as 


-  Sj-i)  +  ijjSj+i  -  Sj) 
Sj+ 1  —  Sj-i 


(116) 


where  j  ^  1  and  N. 

It  can  be  seen  that  KtJ  and  RHSi  are  functions  of  the  coordinates  and 
the  angle  of  attack.  If  we  have  M  airfoil  sections  with  N  control  points, 
the  system  of  equations  (112)  becomes  a  set  of  M  +  N  equations  for  N 
unknowns  7,  and  M  unknowns  V>fc-  M  additional  equations  are  determined 
using  the  Kutta  condition  for  each  airfoil  section.  For  this  purpose  we  use 
an  additional  control  point  at  the  extension  of  the  bisector  of  the  trailing 
edge  and  assume  that  the  streamline  through  other  control  points  of  a 
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given  section  passes  through  this  point,  too.  Applying  eq.  (112)  to  these 
additional  point  Or,  we  can  write  the  Kutta  condition  as  follows 


N 


</>*  +  £  7;  A'ri  =  RHSt 
j= i 


(117) 


Now,  the  complete  system  of  equations  for  an  airfoil  section  can  be  written 
in  matrix  form  as 


'  Ki, i  •  • 

•  KitN  1 

'  7i  ‘ 

■  RHSi  ' 

Kna  •  • 

.  Kn,n  1 

7  N 

— 

RHSn 

Kutta 

condition 

.  . 

.  RHSn+ i  . 

Solution  of  the  system  (118)  gives  us  a  nondimensional  vortex  density  jj 
and  the  stream  function  Now,  the  pressure  coefficient  at  any  point  on 
the  airfoil  surface  may  be  obtained  from  the  velocities  via  the  equation 

Cp(xj,yj)  =  1  -  (rr1-)2  =  1-7 )  (119) 

MX> 

When  this  method  is  extended  to  multi-component  airfoils,  we  have  a  dif¬ 
ferent  value  of  the  stream  function  for  each  airfoil  section.  For  example, 
if  we  have  a  common  configuration  slot  -  main  airfoil  -  flap,  the  matrix 
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equation  can  be  written  as 


K\,x  •  •  Ki,3N  10  0 

7i 

‘  RHSi 

Kn,i  •  •  Kn,3N  10  0 

In 

RHSn 

Kn +i,i  ■  ■  Kn+i,3N  0  10 

7N+1 

RH  Sn+i 

KiN,l  ■  ■  K2N.3N  0  10 

72  N 

— 

RH  S2N 

Kin +1,1  ■  ■  K2N+i,3N  0  0  1 

72N+1 

RHS2N+1 

^3AT,l  •  •  K3N,3N  0  0  1 

7 3N 

RH  Szn 

Kutta  condition  for  slot 

V’i 

RH  S3N+1 

Kutta  condition  for  main  airfoil 

RHS3N+2 

Kutta  condition  for  flap 

.  ^3  . 

RHS3N+3 

k  *  »  «  .ib  '  w  d  i  '-f  j 

(120) 

It  can  be  seen  that  only  right  hand  sides  of  systems  (118)  and  (120)  include 
the  angle  of  attack  a.  If  we  want  to  determine  flow  at  different  angles 
of  attack,  it  is  necessary  to  determine  coefficient  matrix  Kij  once,  and 
recalculate  RHSi  as  a  function  of  the  angle  of  attack. 

9.1  A  DESIGN  METHOD 

In  analysis  mode,  values  of  coordinates  x,  and  y,  and  influence  coefficients 
Kij  are  given  by  the  airfoil  geometry  and  systems  of  equations  (118)  and 
(120)  solved  for  the  surface  velocity  distribution  yj.  However,  for  the  airfoil 
design  problem  either  values  of  surface  velocities  y j  or  pressure  distribution 
are  given  and  the  governing  equations  are  to  be  solved  for  the  airfoil  geom¬ 
etry  Xi  and  y{.  It  is  not  possible  to  obtain  a  direct  solution  for  geometry 
but  instead  we  may  apply  an  iterative  procedure  in  which  the  geometry  of  a 
starting  airfoil  section  is  gradually  modified  until  the  desired  velocity  (pres- 
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sure)  distribution  is  achieved.  Substituting  the  desired  values  for  stream 
function  TpD  and  velocities  in  eqs.  (112)  and  (117)  while  retaining  the 
Kij  from  the  previous  design  iteration  allows  solution  for  the  geometry  of 
a  modified  section.  During  the  iteration  process,  the  x ,•  coordinates  may 
be  kept  constant  and  new  values  of  y,-  calculated.  For  single  airfoil  design, 
the  value  of  a  stream  function  tl>  is  arbitrary  since  its  effect  is  only  to  move 
the  airfoil  up  or  down  relative  to  the  x  —  y  coordinate  system.  However, 
in  the  case  of  multicomponent  airfoils,  the  difference  in  stream  functions 
between  any  two  components  determines  the  flow  through  the  slot  between 
the  components,  and  specification  of  the  stream  function  values  depends  on 
the  type  of  design  (i.e.,  given  geometry  of  one  or  more  sections  -  design  of 
other  components;  or  design  of  all  components).  To  start  the  design  pro¬ 
cess  it  is  necessary  to  define  an  arbitrary  basic  airfoil  shape  to  generate  the 
initial  set  of  influence  coefficients  Kij.  After  that  the  influence  coefficients 
for  each  new  iteration  are  calculated  using  the  coordinates  obtained  from 
previous  iteration.  For  example,  at  iteration  m  of  the  design  procedure  for 
control  points  we  have 


VT  =  —  Mine  +  +  E  V]  i  •  =  1,2— N 

A/  » 


(121) 


and  for  the  trailing  edge  points 


=  77— +  </-?  +  E  KT,  VI 


(122) 


Since  the  location  of  the  trailing  point  is  very  close  to  the  trailing  edge  we 
may  assume  that  y f  is  the  location  of  the  actual  trailing  edge.  In  this  way 
after  each  iteration  we  have  the  location  of  one  point  that  determines  the 
remaining  points  on  the  airfoil  surface.  The  iteration  process  ends  if  either 


i 


r/VtA 

M 


1 


liJ'iViV 


m 

m 

■ 

n 


p©:,: 

$$$ 

'l.l'kVJ 


& 


mm 


It 


y,yAV 

V 

•iVA'A 


condition 


™<ixi[\yr{1  —  yr-1|]  <  (123) 

or 

maic,[|C^-C^-1|]<r!  (124) 

for  tolerances  Ti,r2  is  met.  During  the  design  process,  the  control  points 
of  the  airfoils  are  adjusted  along  vertical  lines  z,  =  const.  This  implies 
that  we  can  not  have  adjustment  of  coordinates  parallel  to  the  freestream 
direction.  In  practical  applications  this  limitation  is  relaxed  because  we 
have  a  prescribed  chord  length  and  we  always  can  have  an  arbitrary  x  - 
coordinate  distribution. 

The  coordinates  y^  that  are  obtained  from  equations  (121  and  122)  are 
panel  control  points.  The  points  on  the  airfoil  surface  can  be  determined 
by  passing  a  curve  through  control  points  and  interpolating  or  projecting  a 
straight  line  through  the  control  points  using  the  trailing  point  as  an  actual 
trailing  edge  of  an  airfoil.  In  the  latter  method,  however,  a  small  error  due 
to  incorrect  location  of  the  trailing  edge  can  propagate  during  iteration  and 
produce  a  saw-shaped  airfoil  that  can  not  satisfy  the  required  pressure  or 
velocity  distribution.  In  order  to  avoid  these  difficulties,  referring  to  fig. 
37a,  after  each  iteration  we  compute  the  actual  position  of  the  trailing  edge 
using  the  first  and  last  control  points  on  the  airfoil,  x-  coordinate  of  the 
trailing  edge  of  the  starting  geometry  and  the  additional  trailing  point.  In 
this  way  we  can  compute  the  coefficient  matrix  Kij  more  accurately  and 
use  the  actual  trailing  edge  as  starting  point  for  determination  of  the  airfoil 
surface.  For  example,  referring  to  fig.  37b,  using  the  equation  of  a  straight 
line  through  two  points  P\  (trailing  edge)  and  C\  (the  first  control  point) 
we  can  determine  location  of  P2  (the  second  point  on  the  airfoil  surface), 
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and  so  on. 

10  DISCUSSION  OF  RESULTS 

10.1  AIRFOIL  ANALYSIS 

For  the  analysis  we  selected  three  characteristic  single  airfoils:  cambered 
Karman  -Treffiz,  cambered  Joukovski  and  a  conventional  NACA  2412  air¬ 
foil.  We  also  consider  multi-component  airfoils. 

•  EXAMPLE  1:  Cambered  Karman  -  Treftz  airfoil  at  a  =  5°  (fig.  38) 

For  this  test  case  we  use  the  stream  function  approach  with  41  panels, 
source  vortex  panel  method  and  higher  order  (  curved  )  panel  method 
with  linear  distribution  of  vorticity  with  60  panels.  It  can  be  seen 
that  the  source  vortex  method  gives  a  larger  pressure  peak  at  the 
leading  edge  on  lower  surface  in  comparison  with  an  exact  solution 
obtained  by  conformal  mapping.In  fact,  neither  surface  singularity 
method  gives  a  completely  accurate  solution  on  the  lower  surface 
,  but  the  curved  panel  method  and  stream  function  approach  can 
handle  leading  and  trailing  edges  approximately  with  the  same  order 
of  accuracy. 

•  EXAMPLE  2:  Cambered  Joukovski  airfoil  at  a  =  2°  (  fig.  39) 

For  analysis  of  this  airfoil  we  use  the  stream  function  approach  with 
41  panels  and  a  higher  order  vortex  panel  method  with  60  panels.  It 
can  be  seen  that  both  methods  give  very  good  agreement  with  exact 
solution.  As  in  the  first  example,  the  cambered  Joukovski  airfoil  gives 
different  velocities  on  upper  and  lower  surfaces  (  high  loading  near  the 
trailing  edge)  that  the  Kutta  condition  (76)  can  not  satisfy.  Neither 
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source  vortex  nor  linear  vortex  panel  method  applied  here  with  this 
Kutta  condition  can  handle  the  present  airfoil  with  given  number  of 
panels,  successfully. 

•  EXAMPLE  3:  NACA  2412  airfoil  at  a  =  8°  (  fig.  40) 

This  test  case  shows  a  comparison  between  linear  vortex  and  stream 
function  panel  methods  with  40  panels.  It  can  be  seen  that  both 
methods  give  approximately  the  same  solution.  The  test  case  also 
proves  that  for  airfoils  with  very  small  camber, trailing  edge  velocities 
are  almost  the  same  and  the  Kutta  condition  (76)  can  be  applied. 

•  EXAMPLE  4:  William’s  configuration  at  a  =  0°  (  fig.  41) 

For  analysis  of  multicomponent  airfoils  the  most  popular  test  case 
is  William’s  configuration  with  flap  at  10°  or  30°.  In  the  present 
example  we  use  the  configuration  with  flap  at  30°.  The  problem  is 
solved  with  41  and  60  panels  on  each  section.lt  can  be  seen  that  a 
better  solution  was  obtained  with  41  panels.  The  solution  on  the  flap 
compares  well  with  the  exact  solution,  while  on  the  lower  surface  of 
the  main  airfoil  there  is  a  small  difference  between  panel  and  exact 
solutions.  The  same  problem  is  observed  by  Seebohm  and  Newman 
[57]  using  a  linear  vortex  panel  method  with  quadratic  extrapolation 
of  the  Kutta  condition.The  present  method  gives  approximately  the 
same  form  of  solution  on  the  main  airfoil  but  gives  a  better  solution 
on  the  flap. 

•  EXAMPLE  5:  Arbitrary  two  component  configuration 

Fig.  42  indicates  an  arbitrary  two  component  configuration  of  two 
NACA  2412  airfoils.  Flap  chord  is  30%  of  the  main  airfoil  chord  with 
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angle  of  deflection  of  10°. The  configuration  is  analysed  at  different 
angles  of  attack:  1°,5°  and  9°.  The  form  of  solution  is  similar  to  that 
obtained  by  Moran,  Cole  and  Wahl  [51]. 


10.2  AIRFOIL  DESIGN 


10.2.1  SINGLE  AIRFOIL  DESIGN 


The  iterative  design  procedure  should  converge  on  an  airfoil  design  which 
gives  exactly  the  required  velocity  distribution.  For  this  purpose  we  apply 
velocity  and  coordinate  convergence  conditions 


Ei  =  |<*  -  uf  |  <  ci 


(125a) 


£2  =  |y,m  -  2/r1 1  <  *2 


(125b) 


where  superscript  m  indicates  iteration  level  and  R  the  required  velocity 
distribution. 


To  start  the  design  process  it  is  necessary  to  specify  a  target  velocity 
distribution.  Convergence  depends  on  how  close  a  starting  geometry  is  to 
the  target  geometry.  If  the  starting  geometry  is  far  from  the  target  design 
process  will  take  a  larger  number  of  iterations.Two  cases  were  tested. 


EXAMPLE  1:  In  this  example  we  attempt  to  design  William’s  main 
airfoil  starting  from  a  NACA  0012  airfoil  (  fig.  43). There  is  a  big 
difference  in  trailing  edges  and  also  a  big  difference  in  pressure  dis¬ 
tribution  in  this  region.  Using  the  known  velocity  distribution,  forty 
panels  were  used  for  discretizing  the  starting  geometry.  After  14  it¬ 
erations  the  velocity  convergence  criterion  (ei  =  0.001)  was  satisfied 
and  the  solutions  differ  only  slightly  at  the  trailing  edge. 
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•  EXAMPLE  2:  NACA  2412  (  fig.  44)  designed  from  a  NACA  0012 
at  the  same  angle  of  attack  (a  =  8°).  The  starting  geometry  was 
very  close  to  the  required  geometry  and  after  9  iterations  the  coordi¬ 
nate  convergence  criterion  was  satisfied.  The  design  process  was  very 
stable  and  very  accurate. 

10.2.2  MULTICOMPONENT  AIRFOIL  DESIGN 

In  multicomponent  airfoil  design  one  can  consider:  (1)  design  of  all  compo¬ 
nents,  (2)  given  main  airfoil  -  design  other  components,  and  (3)  given  flap 
or  slot  -  design  of  the  main  airfoil.  The  second  case  is  the  most  interesting 
and  practical  because  of  several  requirements  that  should  be  satisfied  in 
design  of  high  lift  devices.  These  requirements  are:  (a)  the  flap  and/or  slot 
must  retract  into  the  main  airfoil  to  form  a  good  airfoil  for  cruising  flight, 
( b )  the  deflection  angles  should  provide  reasonable  increase  of  lift,  and  (c) 
the  gap  between  main  section  and  flap  should  provide  enough  energy  to  push 
the  separation  point  of  the  boundary  layer  as  much  as  possible  towards  the 
trailing  edge  of  the  flap. 

•  EXAMPLE  1:  Design  of  William’s  configuration  with  flap  at  30°  ( 
fig.  45) 

The  required  velocity  distribution  was  given  from  the  exact  solution  at 
a:  =  0°.  The  starting  geometry  was  composed  from  a  NACA  0012  as 
main  section  and  from  NACA  0009  as  flap  at  zero  angle  of  deflection. 
It  is  important  that  the  initial  flap  angle  and  spacing  are  not  critical 
because  these  parameters  are  quickly  adjusted  by  the  design  process. 
Since  the  x  -  coordinates  are  not  altered  during  the  design  process, 
the  flap  chord  will  change  if  the  flap  deflection  is  different  from  the 


starting  deflection.  For  example,  if  a  particular  final  flap  chord  is 
desired,  we  can  determine  angle  of  deflection  and  set  starting  flap  at 
that  angle. 

The  required  surface  velocities  on  both  components  were  supplied 
together  with  stream  functions  V’l  and  ^>2  on  the  main  section  and 
flap,  respectively.  Discretization  employed  80  panels  -  40  on  each 
section.  The  design  test  case  is  very  difficult  because  of  a  high  velocity 
peak  at  the  leading  edges  of  both  sections.  Very  small  deviation  in 
coordinates  of  the  leading  edge  will  cause  a  big  change  in  pressure. 
It  can  be  seen  that  pressure  on  the  flap  is  very  close  to  the  required 
pressure  distribution  and  that  at  the  leading  edge  of  the  flap  the 
pressure  peak  is  a  little  overestimated.  On  the  main  section  we  have 
good  agreement  in  pressure  on  the  upper  surface  and  at  the  leading 
and  trailing  edges,  but  underestimated  pressure  on  the  lower  surface. 

EXAMPLE  2:  Here  we  design  the  same  configuration  as  in  the  first 
example  but  using  the  main  William’s  airfoil  and  NACA  0009  as 
starting  geometry  (fig.  46).  That  is,  we  want  to  design  a  flap  with 
a  fixed  main  airfoil.  During  the  design  process  the  velocity  distribu¬ 
tion  on  the  main  airfoil  is  assumed  to  be  unknown  and  is  replaced 
by  the  values  computed  in  the  previous  iteration,  while  geometry  and 
location  of  the  main  section  are  held  constant.  It  can  be  seen  that 
after  five  iterations  the  pressure  distribution  on  the  flap  is  close  to 
the  required  values  but  there  is  a  big  difference  in  pressure  distribu¬ 
tion  on  the  main  section.  The  reason  is  a  relatively  slow  change  of 
pressure  distribution  on  the  flap  that  can  also  be  observed  in  fig.  42 
for  two  NACA  airfoils  at  different  angles  of  attack.  Final  solution 
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was  obtained  after  9  iterations  with  good  agreement  with  pressure 
distribution  in  analysis  mode  (fig.  41) 

11  CONCLUSION 


Surface  singularity  methods  are, a  good  means  for  solution  of  potential  flow 
problems.  They  can  handle  both  single  and  multicomponent  airfoils  in 
analysis  and  design  processes.  The  stream  function  approach  appears  to 
be  very  stable  and  a  simple  panel  method  applicable  for  practical  initial 
analysis  and  design  airfoil  problems. The  method  provides  reduction  in  the 
number  of  panels  in  comparison  with  other  surface  singularity  methods  and 
the  best  results  were  obtained  with  40  -  60  panels  per  airfoil  section. The 
modified  approach  for  computing  the  actual  position  of  the  trailing  edge 
provides  designed  airfoils  with  smooth  shapes  and  stable  design  process. 
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APPENDIX  A 

12  VORTEX  PANEL  METHOD  WITH  LI¬ 
NEAR  VORTEX  STRENGTH 

The  circulation  density  on  each  panel  varies  linearly  from  one  corner  to  the 
other  and  is  continuous  across  the  corner 

7.,-  =  7 j  +  (7j+i  “  7 j)|j  (126) 

Recalling  the  basic  equation  for  the  vortex  panel  method,  the  velocity  po¬ 
tential  at  the  i  -th  control  point  is 

•  y«)  =  VooixiCOsa  -  yisina )  -  ^  /  -^—-tan~1(  — - — )  dsj  (127) 

j=l  *'4j  “  Xcj 

Applying  the  boundary  condition 

=  0  *  =  1,2. ..AT  (128) 

an, 

and  carrying  out  the  differentiation  and  integration  of  equation  (127),  we 
get  a  matrix  equation 

E  =  !•,  (129) 

j+i 

If  we  enforce  the  Kutta  condition  explicitly,  i.e  ,  7i  =  7 n  =  0,  we  get  an 
overdetermined  system  that  can  be  solved  using  the  least-  squares  method. 
The  idea  is  to  find  jj  so  as  to  minimize 

N+l  N 

E2  =  £(Etf.7T>  -&.)2  (130) 

«=1  1=1 

The  minimum  of  (130)  is 

N+l  N 

j-  =  o  =  2E(E^«ii-w, 

t=i  j+i 


(131) 


N  N+l  N+ 1 

13  (  13  KijKimhi  =  13  m  =  1,2. ..AT 

j=i  «=i  t=i 


(132) 


Now,  (132)  represents  a  determined  system.  However,  if  we  use  the  Kutta 
condition  in  the  form  of  eq.  (76),  (127)  can  be  discretized  as  follows 

N 


Z(Klfl1  +  k; f  >)  7i  =  sin(6,  -  a) 

i=i 

Where  the  coefficients  are 

1  =  \bf  +  ag-  kP1 


(133) 


,XK)  „  ,  1  IF  (AD+  BE)G, 

A‘>  -•b  +  5I7 - h - 


(134) 

(135) 


and 


A  =  sin  (0,  -  0j ) ;  B  =  cos  (0,-  —  0y) 

C  =  (xCi  -  X,)2  +  (yc<  -  Ytf 
D  =  (Xj  -  xc,)cos6j  +  (Yj  -  yd)sin6j 
E  =  (xc,  -  X:)sin  0j  -  (yCt  -  Yj)cos  02 
F  =  (n[l  +  i(i?  +  2L,D)] 

G‘At^(crm] 

H  =  (x c,  -  Xj)sin[0i  -  10 2)  +  ( yCi  -  Y2)sin{0i  -  20 2) 
For  i  =  j  the  coefficients  become 
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The  solution  of  the  system  (133)  gives  us  the  unknown  circulation  densities. 
The  tangential  component  of  velocity  can  be  computed  as  follows 


VT,  =  cos  (, 9i  -  a)  +  £(/4Tl)  + 
j= i 


(137) 


where  the  coefficients  for  i  ^  j  are 


and  for  i  =  j 


KlJl]  =  l-AF  -BG-  K&] 

jAt2)  1  HF  ( BD  -  AE)G 
~A+ 2  Lj  +  Lj 


k(Ti)  _  tAT2)  _  £ 


Then,  the  pressure  coefficient  can  be  computed  as 


c„.  =  1  -  vl 


(138) 
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APPENDIX  B 

13  CURVED  PANEL  DISCRETIZATION 


We  use  a  base  coordinate  system  for  an  element:  Referring  to  fig.  (47)  the 
coordinate  system  is  oriented  and  positioned  as  shown. 

a)  £  or  x  coordinate  is  tangent  to  the  curve;  b)  normal  projections 
of  the  x  -  axis  of  the  ends  of  the  curve  S  lie  equal  distance  A  to  the  right 
and  left  (/  =  2A);  c)  tj  -  axis  is  normal  to  the  curve;  d)  general  point  R 
on  the  curve  has  coordinates  (£,77);  and  e)  the  distance  between  some 
arbitrary  point  P(x ,  y )  and  i?(£,  if)  is  r  =  yj(x  —  £)2  -f  (y  —  rj )2 

The  boundary  curve  is  defined  as  function  of  the  £  -  coordinate,  i.e., 
7}  —  viO-  a  neighborhood  of  the  origin  the  curve  and  vorticity  have 
power  series  approximations 


T)  =  of  +  bf  + .  (139) 

-y  _  y°)  _|_  y2)fi2  4- .  (140) 

or  as  a  function  of  panel  curvature 

*  =  z+lc2e  (i4i) 

ds  =  (1  +  2 c2i2)  di  (142) 


Factor  c  depends  on  panel  curvature  and  it  is  obtained  by  three  point  curve 
fitting  using  the  averaged  element  slope 


Cj  = 


-  tan  0j-i)  +  l0^{tan9j+1  -  tanOj) 


/ 1 


n 


(143) 


Using  eq.  (143),  a  second  order  curved  panel  formulation  can  be  written 


—  n,  Vjx 


(144) 


Using  the  Kutta  condition  7X  =  7^,  system  (144)  becomes 

iv-i 

£  tf«7i  = 
i= 2 


(145) 


where  represents  the  left-hand  side  and  RHSi  the  right-hand  side  of 
system  (144).  The  influence  coefficients  for  normal  velocity  component  can 
be  expressed  as 


tAN)  _  j-ANr)  k(n2) 


(146) 


<■’  -  4’ <5  +  f  +  ffad  +2c?f7£MC  (147) 


“  /-4  (5  +  +  +  <148> 


Where 

50  (&  -  Y,  )  sin  0,-  +  (s(Q  -  -Xj)  cos  0,- 

(5nj,J  (x,  -  X,-)*  +  (tt V,)2  U  j 

Xj  =  Xj  +  £cos  6j  —  Cj^sin  6j 
Yj  =  yj  -f  (sin  0j  —  Cjfcos  0; 

Coordinates  Xj,Yj  are  new  control  points  obtained  by  parabolic  curve  ap¬ 
proximation  for  each  panel,  and 

Kij  =  1  for  i  =  j 
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Kij  =  -1  for  i  =  j  +  1 

The  influence  coefficient  for  tangential  velocity  component  can  be  written 


k!P  =  k\Jx)  +  k\J2) 


■c3)(i  +  2cJ2_1e2)(£),J-irfc 


C3/3 

f°r  i  =  j 


c3  /3 

cj-ib-x 


/or  i  =  j  +  1 


(150) 


=  JP <r +  £  ■ +  (151) 

— (y-  _  Yj)  cos  6i  +  ( X(i )  -  Xj )  sin  6i 


Xj)  +  (y,  -  Yjf 
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NACA  0012 


ELIPTIC  SCHEME 


xx  +  ^  xx  “  ^ 

^yy  +  1lyy  =Q(^1) 


GRID  SIZE 

ERROR 

ITER 

27  x  20 

0.0001 

78 

27  x20 

0.0001 

71 

27  x  20 

0.0001 

68 

27  x20 

0.0001 

57 

61  x  28 

0.0001 

99 

61  x28 

0.001 

479 

TWO  NACA  0012  AIRFOILS  -  FLAP  ANGLE  25° 


69x20 

0.0001 

275 

69x20 

0.001 

116 

CPU 

time  (sec) 


11.722 

10.708 

10.264 

8.649 

48.561 

239.259 


NACA  0012  **  PARABOLIC  SCHEME  ONLY  ( R  =  3) 

GRID  SIZE 

N°  OF  GRID  POINTS 

CPU  time 

61  x  28 

1708 

0.252 

80x25 

2000 

0.288 

80x40 

3200 

0.469 

100  x  40 

4000 

0.571 

120x80 

9600 

1.354 

120  x  100 

t 

12000 

1.686 

PARABOLIC 

SCHEME 


CPU 

time  (sec) 


0.083 

0.088 

0.090 

0.090 

0.288 

0.283 


1.80 

105.702 

0.340 

1.81 

48.301 

0.342 

Table  1  -  Comparision  of  the  parabolic  and  elliptic  schemes 
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Fig.  1  Discretization  methods 
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Fig.  2  Grid  generation  schemes 


Fig.  5  Grid  spacing  control:  a)  The  effect  of  the  Laplacian; 

b)  the  effect  of  the  forcing  functions  P  and  Q; 

c)  the  effect  of  the  amplitude  factors  in  the  forcing 
functions 
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Fig.  10  C-grid  about  an  NACA  0012  airfoil  generated 
by  the  Laplace  system 
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Fig.21  O-grid  about  an  NACA  0012  airfoil  generated  by  the 
Parabolic  scheme  with  strong  clustering  at  the  airfoil 
surface 

(grid  size  41  x  20,  CPU  time  0.207  sec.) 
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Fig.  23  H-grid  generated  by  the  parabolic  scheme 
with  strong  clustering  and  orthogonality 
at  the  airfoil  surface 
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Fig.  24  C-grid  about  an  NACA  0012  airfoil 
generated  by  the  parabolic  scheme 


m 


iVjViV 


l*KWWi 

p& 

wM 

KvhWv 
W-‘Am 


fiS® 

kW'W 

tv.'iW 

wivivi* 
stem 


wm 

slviv!' 

fjtM- 

rlvlvlvi 

if  *HV 


Fig.  36  Geometry  definition  for  influence  coefficients  K  jj 


Fig.  37  a)  Determination  of  the  actual  position  of  the 
trailing  edge; 

b)  Determination  of  surface  points 
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Fig.  46  Design  of  William's  configuration  from  an  NACA  0012 
airfoil  and  fixed  main  airfoil  geometry 
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